F49620-03- 1-0204 


Final  Report 


NON-SYNCHRONOUS  VIBRATION  OF  TURBOMACHINERY  AIRFOILS 

AFOSR  Contract  F49620-03- 1-0204 


DUKE 


EDMUND  T.  PRATT,  JR. 

SCHOOL  OF 
ENGINEERING 


Final  Report 

By 

Robert  E.  Kielb,  Senior  Research  Scientist  and  Principal  Investigator 
Kenneth  C.  Hall,  Professor  &  Chair 
Meredith  Spiker,  Graduate  Student 
Jeffrey  P.  Thomas,  Research  Assistant  Professor 

Department  of  Mechanical  Engineering  and  Materials  Science 
Edmund  T.  Pratt,  Jr.  School  of  Engineering 
Duke  University 
Durham,  NC  27708-0300 


Lt.  Col.  Rhett  Jeffries 
AFOSR  Program  Manager 
AFOSR/NA 

875  North  Randolph  Street 
Suite  325,  Room  3112 
Arlington,  VA  22203-1768 


Duke  University 


Page  1 


March  2006 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection 
of  information,  including  suggestions  for  reducing  the  burden,  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports 
(0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be 
subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY) 

30-05-2006 


4.  TITLE  AND  SUBTITLE 


2.  REPORT  TYPE 


Final  Performance 


Non-Synchronous  Vibration  of  Turbomachinery  Airfoils 


3.  DATES  COVERED  t From  -  To) 

01-04-2003  to  28-02-2006 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 

F49620-03- 1-0204 


5c.  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S) 

Kielb,  Robert,  E.  Senior  Research  Scientist  and  Principal  Investigator 
Hall,  Kenneth,  C,  Professor  &  Chair 
Spiker,  Meredith,  A,  Graduate  Student 
Thomas,  Jeffrey,  P.  Research  Assistant  Professor 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Department  of  Mechanical  Engineering  and  Materials  Science 
Edmund  T.  Pratt,  Jr.  School  of  Engineering 
Duke  University 
Durham,  NC  27008-0300 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Lt.  Col.  Rhett  Jeffries 
AFOSR  Program  Manager 
875  Noith  Randolph  Street 
Suite  325,  Room  3112 
Arlington,  VA  22203-1768 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 


11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

AFRL-SR-AR-TR-06-0269 


14.  ABSTRACT 

The  goal  of  this  research  was  to  develop  a  fast  computational  method  to  be  able  to  understand  non-synchronous  vibrations  in 
turbomachinery  blades  in  the  design  stage.  The  numerical  approach  involved  the  use  of  a  two-  and  three-dimensional  nonlinear 
Navier- 

Stokes  unsteady  harmonic  balance  method.  As  an  initial  demonstration  of  the  method,  the  flow  over  a  cylinder  and  a  two- 
dimensional 

airfoil,  the  2D  Cl  case,  were  investigated  and  results  showed  good  agreement.  For  the  3D  Cl  case  and  a  running  tip  clearance  of 
0.020  in., 

two  irrational  frequencies  were  present  in  the  solution  so  the  current  HB  method  could  not  be  used.  In  order  to  justify  the  merits  of 
the  HB 


15.  SUBJECT  TERMS 

turbomachinery,  blade  vibrations,  flutter,  harmonic  balance 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

ABSTRACT 

OF 

PAGES 

Uncalssified 

Unclassified 

Unclassified 

72 

Robert  Kielb 

19b.  TELEPHONE  NUMBER  (Include  area  code) 
919-660-5327 


Standard  Form  298  (Rev.  8/98) 
Prescribed  by  ANSI  Std.  Z39.18 


F49620-03- 1-0204 


Final  Report 


TABLE  OF  CONTENTS 

Page 

1.  INTRODUCTION . 4 

2.  THEORY  . 11 

Governing  Equations . 13 

Non-Dimensionalization . 16 

Application  of  Harmonic  Balance  Method  to  Governing  Equations . 17 

Utilization  of  Time  Domain  Methodology . 19 

Numerical  Solution  Technique . 19 

Grid  Generation . 20 

Lax-Wendroff  Method . 20 

Application  of  Initial  and  Boundary  Conditions . 22 

Frequency  Search  Procedure . 22 

Phase  Error  Method . 23 

3.  APPLICATIONS  OF  THE  HARMONIC  BALANCE  METHOD . 24 

3.1  Stationary  Cylinder  in  Cross  Flow . 24 

Application  of  Phase  Error  Method . 25 

Strouhal  -  Reynolds  Number  Relationship . 26 

Stationary  Cylinder  Flow  Simulation . 28 

Unsteady  Cylinder  Lift . 30 

Experimental  Procedure  and  Results . 32 

Mean  Cylinder  Drag . 33 

3.2  Cylinder  with  Prescribed  Motion . 34 


Duke  University 


Page  2 


March  2006 


F49620-03- 1-0204 


Final  Report 


Technique  Used  to  Determine  Bounds  of  Lock-In  Region . 35 

Effect  of  Mesh  Refinement . 37 

Addition  of  Filter  to  Reduce  Far  Field  Boundary  Effects . 38 

Unsteady  Lift  for  Prescribed  Motion . 39 

Real  and  Imaginary  Part  of  the  Unsteady  Lift . 39 

Amplitude  of  the  Unsteady  Lift . 41 

Phase  Shift  as  a  Function  of  Strouhal  Number . 43 

Aerodynamic  Damping . 44 

Mean  Cylinder  Drag  for  Prescribed  Motion . 45 

3.3  Elastically  Mounted  Cylinder . 46 

Aeroelastic  Cylinder  Model . 47 

Non-Dimensional  Parameters . 48 

Application  of  Newton-Raphson  Technique . 50 

Experimental  Procedure  and  Results . 52 

3.4  Two-Dimensional  Airfoil  (Cl  Case) . 53 

3.5  Three  Dimensional  Cases . 57 

Cl  Case . 57 

HI  Case . 67 

H2  Case . 68 

4.  SUMMARY  AND  CONCLUSIONS . 68 

REFERENCES . 71 


Duke  University 


Page  3 


March  2006 


F49620-03- 1-0204 


Final  Report 


1.  INTRODUCTION 

Currently,  in  turbomachinery,  the  two  main  design  considerations  are  blade  forced 
response  and  flutter.  Forced  response  of  the  blades  result  due  to  aerodynamic  excitations 
from  upstream  wakes  and  occur  at  integer  multiples  of  the  vane  passing  frequency. 

Flutter  is  a  self-excited  aeroelastic  instability  in  which  the  aerodynamic  forces  on  the 
blade  add  to  the  blade  vibration  itself.  In  both  these  cases,  the  vibrations  can  be  large 
enough  to  cause  the  blades  to  fracture.  However,  in  recent  years,  engine  manufacturers 
have  encountered  a  new  class  of  aeromechanical  problems.  This  category  generally 
includes  separated  flow  vibrations  (SFV)  and  non-synchronous  vibrations  (NSV). 
Separated  flow  vibrations  are  similar  to  buffeting  in  wings  and  occur  when  an  unsteady, 
separated  flow  is  generated  over  a  row  of  turbine  blades  that  randomly  excites  blade 
vibration  (Hall  2004,  675).  In  addition,  they  generally  have  a  broadband  frequency 
content.  Non-synchronous  vibrations  are  similar  to  separated  flow  vibrations  except  they 
are  usually  well-ordered,  can  occur  away  from  a  stalled  condition,  and  occur  at  one 
dominant  frequency.  The  blade  vibrations  are  generally  frequency  and  phase  locked. 
Typical  non-turbulent  examples  of  these  phenomena  include  the  vortex  shedding  flow 
over  a  car  antenna,  power  lines,  cables,  or  off-shore  risers  and  all  can  cause  large 
amplitude  vibrations.  Although  both  non-synchronous  vibrations  and  separated  flow 
vibrations  are  emerging  research  areas,  this  effort  is  only  concerned  with  the  study  and 
prevention  of  non-synchronous  vibrations. 

Non-synchronous  vibrations  in  turbine  engine  blades  are  the  result  of  the 
interaction  of  an  aerodynamic  instability,  such  as  vortex  shedding,  with  the  vibrations  of 
the  blades.  It  is  a  serious  aeroelastic  problem  that  has  been  observed  by  most  engine 
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companies  and  can  ultimately  lead  to  blade  fracture.  Currently,  there  is  limited 
knowledge  of  this  phenomenon,  so  failures  typically  do  not  occur  until  the  testing  phase. 
As  a  result,  engineers  must  redesign  the  blade  or  the  rotor.  This  procedure  is  costly  and 
greatly  increases  production  time. 

Therefore,  the  overall  goal  of  the  NSV  research  effort  at  Duke  is  to  develop  a  fast, 
systematic  method  for  engine  companies  to  prevent  and  eliminate  NSV.  This  will  be 
accomplished  by  using  a  novel  computational  fluid  dynamic  technique,  namely  the 
nonlinear,  unsteady  harmonic  balance  method  developed  by  Hall,  et  al.  to  solve  the  three- 
dimensional  Reynolds  averaged  Navier-Stokes  equations  governing  the  flow  through  an 
engine  blade  passage  (Hall,  Thomas,  and  Clark  2002,  879).  This  method  offers  distinct 
computational  advantages  over  currently  used  methods.  In  particular,  it  requires  one  to 
two  orders  of  magnitude  less  computational  time  than  conventional  time  marching 
computational  fluid  dynamic  (CFD)  methods  (Hall,  Thomas,  and  Clark  2002,  886). 
Furthermore,  the  harmonic  balance  method  can  handle  large  flow  disturbances  while 
time-linearized  equations  are  only  valid  for  small  unsteady  perturbations.  However, 
typically  in  a  flow  analysis,  the  oscillation  frequency  is  known  beforehand  and  therefore 
serves  as  an  input  to  the  system.  Since  the  goal  of  the  research  is  to  find  the  NSV 
frequency,  a  new  phase  error  method  is  employed  to  achieve  quick  solution  convergence. 
This  novel  frequency  search  technique  utilizes  the  phase  difference  between  successive 
iterations  of  an  unsteady  first  harmonic  quantity  such  as  unsteady  lift  and  was  developed 
by  Kenneth  Hall  (private  communication).  It  requires  only  a  few  harmonic  balance 
solution  computations  to  determine  the  precise  NSV  frequency.  This  method  will  be 
described  in  detail  in  the  latter  portion  of  the  theory  section.  Therefore,  by  integrating 
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these  numerical  solution  techniques,  a  more  computationally  efficient  method  for 
predicting  both  the  NSV  frequency  and  amplitude  will  be  developed. 

An  accurate  prediction  of  the  NSV  is  important  because  it  provides  useful 
Campbell  diagram  information  for  rotating  blade  design  purposes.  A  hypothetical 
Campbell  diagram  is  presented  in  Figure  1.  The  Campbell  diagram  is  a  plot  of  the 
blade’s  natural  frequencies  of  vibration  in  Hz  as  a  function  of  the  rotational  speed  of  the 
engine  in  RPM.  The  diagonal  lines  of  constant  slope  represent  engine  order  excitations 
lines,  i.e.  multiples  of  the  operation  speed  or  machine  harmonics.  The  nearly  horizontal 
lines  represent  natural  frequencies  of  particular  blade  modes,  i.e.  first  flex  (IF)  or  first 
torsion  (IT).  As  a  result  of  the  harmonic  balance  analysis,  the  NSV  frequencies  could  be 
added  to  the  Campbell  diagram  to  identify  regions  where  the  frequency  of  the  flow 
instability  might  interact  with  a  blade  natural  frequency,  resulting  in  blade  vibrations. 
These  possible  points  of  intersection  have  been  identified  on  the  Campbell  diagram. 
Presently,  engineers  commonly  use  the  Campbell  diagram  to  detennine  potential  resonant 
response  problems  by  examining  the  crossing  synchronous  excitations  (such  as  vane 
passings)  with  the  blade  natural  frequencies.  In  addition,  the  presence  of  asynchronous 
vibrations  such  as  flutter  can  be  determined  by  examining  the  Campbell  diagram.  For  all 
cases,  the  stress  levels  at  the  crossing  points  can  be  high  enough  to  cause  high  cycle 
fatigue  of  the  blades  and  therefore,  greatly  reduce  the  component’s  life  span.  However, 
other  aeromechanic al  problems  such  as  SFV  and  NSV  are  not  presently  placed  on  a 
Campbell  diagram  because  they  are  not  well  enough  understood  and  a  design  analysis 
method  has  not  yet  been  developed.  Therefore,  the  proposed  method,  when  extended  to 
blades,  will  allow  engineers  to  better  understand  NSV  behavior  and  to  predict  the 
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occurrence  of  vortex-induced  vibrations.  Once  this  fluid  dynamic  instability  becomes 
better  understood,  it  will  be  possible  to  couple  the  NSV  with  blade  motion. 


Figure  1.  Hypothetical  Campbell  Diagram  Demonstrating  Design  Considerations 
Including  Non-Synchronous  Vibrations  (NSV) 

As  an  initial  step  towards  better  understanding  the  NSV  in  turbine  engine 
configurations,  it  is  advantageous  to  investigate  a  well-known  two-dimensional  case  that 
exhibits  similar  NSV  features  such  as  shedding  flow  about  a  circular  cylinder.  In  this 
case,  the  flow  separates  from  the  cylinder’s  surface  and  results  in  a  wake  region 
characterized  by  a  low  pressure  behind  the  body  (Fox  and  McDonald  1992,  34).  If  a 
cylinder  is  placed  in  a  low  Reynolds  number  flow  (Re>47),  vortices  are  shed  alternately 
from  the  top  and  bottom  of  the  cylinder.  This  phenomenon  is  generally  known  as  Von 
Karman  vortex  shedding. 

The  flow  behind  the  cylinder  becomes  three-dimensional  for  Reynolds  numbers 
greater  than  180.  This  case  is  much  more  difficult  to  model,  so  this  preliminary  study 
only  considers  flow  in  the  Reynolds  number  range  of  47  <  Re  <  180.  The  Reynolds 
number  measures  the  relative  significance  of  the  viscous  effects  as  compared  to  inertial 
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effects.  At  very  low  Reynolds  number  (i.e.  Re  <  5),  the  flow  remains  attached  to  the 
cylinder.  As  the  Reynolds  number  is  increased,  the  flow  starts  to  separate  from  the 
surface  of  the  cylinder  (5  <  Re  <  47).  At  a  Reynolds  number  greater  than  47,  the  flow 
begins  to  shed.  If  the  Re  is  continually  increased,  the  vortices  begin  to  detach  themselves 
from  the  cylinder,  resulting  in  a  periodic  wake  with  two  staggered  rows  of  vortices.  This 
regular  pattern  of  vortices  is  known  as  a  Von  Karman  vortex  street  (see  Figure  2). 


Figure  2.  Development  of  Von  Karman  Vortex  Streets  in  the  Laminar  Flow  Regime 

for  a  Re  =  150  (Williamson  2004,  480) 

The  flow  over  a  stationary  cylinder  can  be  determined  by  two  non-dimensional 
parameters  -  the  Reynolds  number  and  the  Strouhal  number.  For  a  given  Reynolds 
number,  the  flow  sheds  at  a  distinct  frequency.  The  Reynolds  number  represents  a  non- 
dimensional  ratio  of  the  inertial  forces  in  a  fluid  to  the  viscous  forces.  It  is  given  by: 


Re  = 


UD 


(1) 


where  Uoo  is  the  velocity  of  the  fluid,  D  is  the  characteristic  dimension  (i.e.  diameter  of 
the  cylinder),  and  Voo  is  the  kinematic  viscosity  of  the  fluid.  Whether  the  flow  is  laminar 
or  turbulent  is  primarily  governed  by  the  Reynolds  number.  This  study  is  concerned  with 
the  laminar  flow  over  a  circular  cylinder.  Another  important  parameter  related  to  the 
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flow  over  a  cylinder  is  the  Strouhal  number.  This  is  the  non-dimensional  vortex 
shedding  frequency.  It  is: 

St  =  (2) 

U 

oo 

where  fst  is  the  dimensional  vortex  shedding  frequency,  D  is  the  cylinder  diameter,  and 
Uoo  is  the  steady  velocity  of  the  uniform  flow.  A  relationship  between  Reynolds  number 
and  Strouhal  number  was  determined. 

Another  important  measure  of  vortex  shedding  is  the  unsteady  lift  on  a  stationary 
cylinder.  Lift  is  primarily  due  to  variations  in  pressure  on  the  cylinder’s  surface  and  is 
the  integrated  result  of  the  pressure  loading  on  the  cylinder  (Norberg  2001,  460).  Based 
on  experimental  results  and  2-D  simulations,  it  can  be  seen  that  the  amplitude  of  the  RMS 
lift  coefficient  increases  rapidly  within  the  laminar  shedding  regime  from  the  onset  of 
shedding  at  a  Reynolds  number  of  47  (see  Figure  10).  The  onset  can  be  characterized  as 
a  supercritical  Hopf  bifurcation  and  is  typically  described  by  the  Stuart-Landau  equation 
(Norberg  2001,  464).  The  unsteady  lift  can  be  used  to  determine  the  transition  from 
steady  flow  to  Von  Kannan  vortex  shedding  in  the  two-dimensional  wake  behind  a 
circular  cylinder  in  low  Reynolds  number  flow. 

Finally,  the  HB  methodology  was  applied  to  both  two-dimensional  and  three- 
dimensional  real  world  turbomachinery  blade  applications.  Numerous  other  researchers 
have  studied  the  NSV  problem  in  turbomachinery  and  one  particular  area  that  has 
received  considerable  interest  is  tip  flow  instability.  Examples  of  work  in  this  area  can  be 
found  in  Mailach  (1999),  Mailach  et  al.  (2000  &  2001),  Marz  et  al.  (1999  &  2001),  Inoue 
et  al.  (1999),  Lenglin  &  Tan  (2002),  and  Vo  (2001).  Mailach  et  al.  (2001)  present  results 
from  both  a  four  stage  low  speed  research  compressor  and  a  linear  cascade,  and  conclude 
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that  the  tip  flow  instability  is  a  vortex  interaction  effect  that  produces  a  multi-cell 
circumferentially  traveling  wave.  This  phenomenon  is  found  near  the  stall  line  with  a 
relatively  large  tip  clearance  (greater  than  2%  of  tip  chord).  A  new  Strouhal-type  number 
is  proposed  to  characterize  the  frequency  of  the  oscillation.  Marz  et  al.  (2001)  also  found 
a  tip  flow  instability  on  a  low  speed  fan  rig  near  the  stall  line  and  with  a  large  tip 
clearance.  This  paper  also  presents  the  results  of  an  unsteady  CFD  model  that  predicts  a 
frequency  of  950  Hz  compared  to  the  measured  value  of  880  Hz. 

Camp  (1999)  reports  that  this  vexing  problem  has  been  observed  in  a  high  speed 
compressor.  As  a  result,  an  experimental  study  was  performed  in  a  low  speed 
compressor  facility.  A  helical  acoustic  structure  (circumferential  Mach  number  of  0.84) 
was  found  using  casing  dynamic  pressure  transducers.  It  was  also  found  that  there  were 
step  changes  in  response  frequency  as  the  flow  rate  was  changed.  Although  not  proven, 
the  author  speculates  that  the  phenomenon  involves  vortex  shedding  from  the  blades  that 
excites  the  helical  acoustic  cavity  modes,  which,  in  turn,  excite  the  blades. 

In  this  study,  we  examined  three  different  engine  configurations.  Based  on  input 
from  industry,  NSV  was  encountered  in  experimental  rig  testing  for  two  of  the  three  test 
cases.  In  particular,  we  studied  a  modem  first  stage  compressor  rotor  blade  (Cl),  a 
modem  first  stage  fan  blade  (HI),  and  a  modem  fan  vane  blade  (H2).  Although  NSV 
was  not  encountered  for  the  HI  case,  an  analysis  was  performed  to  show  that  NSV  is  not 
predicted.  As  an  initial  step,  a  flow  instability  about  a  two-dimensional  airfoil  section  of 
the  Cl  blade  was  examined  at  an  off-design  condition.  Finally,  a  full  three-dimensional 
analysis  of  all  three  blades  was  performed  to  predict  both  the  frequency  and  amplitude  of 
the  fluid  dynamic  instability. 
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2.  THEORY 

In  general,  the  unsteady  perturbations  are  significant  and  serve  as  a  useful  test 
case  for  the  verification  of  the  harmonic  balance  method  as  a  valuable  tool  for  solving 
NSV  problems  in  turbomachinery.  The  flow  over  the  cylinder  is  modeled  by  the  two- 
dimensional  nonlinear,  unsteady  Navier-Stokes  equations.  Two  currently  available 
analysis  techniques,  namely  the  time  domain  approach  and  the  time-linearized  frequency 
domain  approach,  were  considered  and  the  merits  of  each  were  examined.  Finally,  the 
hannonic  balance  method  developed  by  Kenneth  Hall  will  be  described.  For  time- 
periodic  flows,  Hall’s  method  can  be  computationally  more  efficient  than  the  other 
approaches  and  it  is  not  limited  by  the  assumptions  required  by  a  linearized  solution. 

In  the  past,  most  researchers  have  employed  a  time-domain  approach  to  model  the 
fluid  flow  over  bluff  bodies  (Ni  1987,  He  1993).  In  this  case,  the  governing  Navier- 
Stokes  equations  are  discretized  over  the  computational  grid  and  then  the  flow  solution  is 
marched  in  time  using  conventional  CFD  techniques  and  applying  appropriate  unsteady 
boundary  conditions.  The  discretization  of  the  flow  equations  requires  the  solution  of  a 
set  of  coupled  nonlinear  equations  at  each  physical  time  step.  Commonly,  an  inner 
iteration  is  introduced  with  a  pseudotime  variable  and  marched  to  a  steady  state.  The 
converged  solutions  obtained  at  the  end  of  the  inner  iteration  represents  a  solution  of  the 
equations  at  the  end  of  the  physical  time  step  (McMullen,  Jameson,  and  Alonso  2001,  1). 
Therefore,  this  procedure  must  be  repeated  numerous  times  to  achieve  one  physical  time 
step  and  resolve  the  features  of  one  period  in  the  oscillation  of  the  flow.  Many 
convergence  acceleration  techniques  such  as  multi-grid  and  variable  pseudotime  stepping 
can  be  implemented  in  the  inner  iteration.  However,  as  more  periods  are  needed  for 
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convergence,  the  computations  can  become  extremely  costly.  Although  this  approach  is 
straightforward  to  implement  and  can  model  both  nonlinear  as  well  as  linear 
disturbances,  it  typically  takes  considerable  computational  time  to  ensure  numerical 
stability.  Therefore,  this  approach  is  too  expensive  for  routine  use  by  industry  to  predict 
and  design  for  NSV. 

Another  popular  technique  for  modeling  the  flow  over  turbomachinery  blades  is 
the  time-linearized  frequency  domain  approach.  This  method  has  been  used  to  model 
transonic  flow  problems  as  well  as  two-  and  three-dimensional  viscous  flows  in 
turbomachinery  (Clark  and  Hall  2000,  Hall  and  Crawley  1989).  In  this  method,  the 
problem  is  divided  into  two  parts:  the  steady  flow  and  the  unsteady  perturbations  in  the 
flow.  The  steady,  or  time-mean,  part  of  the  flow  is  solved  using  well-known  CFD 
techniques  and  the  unsteady  part  is  assumed  to  be  small  and  harmonic  in  time  (cK°l).  The 
governing  Navier-Stokes  equations  can  then  be  linearized  about  the  mean  flow  solution 
to  determine  a  set  of  equations  that  describe  the  unsteady  flow  component  and  they  can 
be  solved  relatively  easily.  By  replacing  the  time  derivative  by  jco  where  co  is  the 
frequency  of  the  unsteady  flow  disturbance,  the  resulting  time-linearized  equations  can  be 
solved  inexpensively.  A  solution  to  the  complex  hannonic  amplitude  is  obtained  at  a 
given  frequency.  When  compared  to  unsteady  time  marching  methods,  the  time- 
linearized  frequency  domain  method  is  much  more  efficient  but  it  is  limited  by  the  linear 
assumption  of  unsteadiness. 

This  study  utilizes  the  harmonic  balance  technique  for  nonlinear,  unsteady  flows 
developed  at  Duke  University  by  Kenneth  Hall  in  2002.  It  was  originally  developed  as  a 
mathematical  tool  to  study  the  behavior  of  harmonic  ordinary  differential  equations 
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(Maple  2002,  3).  However,  in  recent  years,  Hall  extended  the  application  to  the  study  of 
turbomachinery  flows.  In  this  novel  method,  the  unsteady  flow  is  assumed  to  be  periodic 
time  and  can  be  represented  by  a  Fourier  series  in  time  with  frequencies  that  are  integer 
multiples  of  the  fundamental  excitation  frequency.  Then,  using  the  hannonic  balance 
technique,  a  set  of  coupled  partial  differential  equations  is  obtained  for  the  unknown 
Fourier  coefficients  of  the  flowfield.  A  source  term,  dU/dt,  accounts  for  the  unsteadiness. 
Finally,  a  pseudo-time  tenn  is  introduced  into  the  hannonic  balance  equations  and  the 
solution  is  marched  like  a  steady  CFD  solver.  From  this  result,  the  lift,  drag,  etc.  for  each 
harmonic  can  be  determined  from  the  calculated  Fourier  coefficients.  For  a  given 
problem,  it  is  desirable  to  use  the  minimum  number  of  hannonics  because  the 
computational  cost  is  proportional  to  the  number  of  harmonics  included.  Local  time 
stepping,  preconditioning,  and  multi-grid  can  be  use  to  accelerate  convergence. 

Therefore,  a  key  aspect  of  the  methodology  is  its  ability  to  exploit  the  use  of  conventional 
CFD  techniques.  The  theory  presented  in  Hall’s  paper  entitled,  “Computation  of 
Unsteady  Flow  in  Cascades  Using  a  Hannonic  Balance  Technique”,  is  summarized  in  the 
following  sections. 

Governing  Equations 

The  flow  over  a  cylinder  can  be  best  represented  by  the  Reynolds-averaged,  two- 
dimensional  Navier-Stokes  equations  for  unsteady,  nonlinear  flow.  The  derivation  of 
these  equations  is  well  documented  in  literature  and  arises  from  applying  three 
conservation  laws:  conservation  of  mass,  conservation  of  momentum,  and  conservation 
of  energy  (Anderson  1984,  181).  The  flow  is  assumed  to  be  compressible  and  viscous. 
Therefore,  the  resulting  equations  of  motion  in  integral  fonn  are: 
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d_ 

dt 


JJ  Udxdy 


+ 


(F  -  U  —)dy  - 

'sdf  dt 


go  ( G-U ^)dx  =  |£  Sdxdy 


(3) 


where  Dc  is  a  deforming  control  volume  bounded  by  the  control  surface  dDc.  U  is  a 
vector  of  conservation  variables,  F  and  G  are  flux  vectors  in  the  x  and  y  directions,  and  S 
is  the  source  vector.  The  quantities  dfldt  and  dg/dt  are  the  x  and  y  components  of  the 
velocity  of  the  control  surface  dDc.  Based  on  the  strong  conservative  law  form  of  the 
Navier-Stokes  equations  given  above,  the  vector  U  and  the  inviscid  and  viscous 
components  of  F,  G  and  S  can  be  written  in  vector  fonn  as: 
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By  examining  the  governing  equations  in  compact  vector  form,  the  first  row  corresponds 
to  the  continuity  equation,  the  second  and  third  rows  are  the  momentum  equation  in  the  x 
and  y  directions  respectively,  and  the  final  equation  is  the  conservation  of  energy 
equation.  In  these  equations,  p  is  the  density  of  the  fluid;  u  and  v  are  the  velocity 
components  in  the  x  and  y  directions,  respectively;  e  is  the  total  internal  energy;  h  is  the 
total  enthalpy;  and  p  is  the  static  pressure.  If  the  divergence  theorem  is  applied,  a 
differential  form  of  the  Navier-Stokes  is  obtained: 

8U  |  dF(U)  |  dG(U)_s 

dt  dx  dy 
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Assuming  the  fluid  (air  in  this  case)  is  an  ideal  gas  with  constant  specific  heats,  it  is 
possible  to  express  the  pressure  and  enthalpy  in  terms  of  the  conservation  variables.  The 
equation  of  state  is  given  by: 

P  =  pRT  (6) 

where  T  is  the  temperature  and  R  is  the  universal  gas  constant.  In  addition,  for  a  perfect 
gas,  the  following  relationships  exist: 


e  =  cvT ,  h  =  c  T ,  /  =  —  ,  cv  = 


R 


/-I 


S  = 


jR_ 

7-1 


(V) 


where  cv  is  the  specific  heat  at  constant  volume,  cp  is  the  specific  heat  at  constant 
pressure,  and  7  is  defined  as  the  ratio  of  these  two  specific  heats.  Therefore,  the  pressure 
and  enthalpy  can  be  determined  and  are  given  by: 

P  =  (7  - 1)  {pe ~  1  /  2  p[{puf  +  (pv) 2  ] }  (8) 


h  =  ( pe  +  p)/  p 


(9) 


The  shear  stresses  txx,  rxy,  tvv  can  be  written  in  terms  of  the  viscosity  and  the  shear  rate  in 
the  x  and  y  directions,  respectively: 


=  fA 


f  \  du  2  dv  ' 
3  dx  3  dy  y 


x  =  u\ 

xv  r*\ 


du  dv 

- 1 - 

dry  dx 


Tyy  =  M 


f  4  dv  2  du  N 

3  dy  3  dx  y 


(10) 


(ID 


(12) 


where  //  is  the  molecular  viscosity.  Furthermore,  the  tenns  ihx  and  Thy  in  the  energy 
equation  are  given  by: 

=UT*X+VTXy  “9,  (13) 
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r,  =ut  +  vr  —a 

hy  xy  yy  *  y 

where  qx  and  qy  are  the  x  and  y  components  of  the  heat  flux  and  can  be  written  as: 


ycp  )dT 


Prz  J  dx 


VCP  |  dT 

P  h  j  fy 


where  Prp  is  the  laminar  Prandtl  number  and  ffT/dx  and  dT/dy  are  the  temperature 
gradients  in  the  x  and  y  directions,  respectively.  The  variable  //  is  detennined  from 
Sutherland’s  Law  and  results  from  the  use  of  a  dynamic  molecular  viscosity. 
Sutherland’s  formula  for  viscosity  is  given  by: 


/  rji  /  rri  \3/2  ^  S0 

m  =  Mt/t0)  — — — 

\*0  +  *^0 


where  jio  =  3.583  x  10'7  at  a  reference  temperature  T0  =  491.4°R  and  the  constant  So  = 
199.8°R.  For  air  at  standard  conditions,  PrL  =  0.72.  Therefore,  by  using  knowledge  of 
the  flow  and  the  principles  of  thermodynamics,  it  is  now  possible  to  implement  the 
hannonic  balance  method  to  solve  the  governing  two-dimensional,  nonlinear,  unsteady 
Navier-Stokes  equations. 

Non-Dimensionalization 

In  order  to  make  the  computations  easier,  the  equations  are  often  put  into  a  non- 
dimensional  form.  This  allows  the  system  parameters  such  as  Reynolds  number,  Mach 
number,  etc.  to  be  varied  independently  (Anderson  1984,  191).  Therefore,  the  variables 
of  interest  were  non-dimensionalized  as  follows: 


,  *  ,  y  ,  u 

x  = - ,  V  = - ,  u  = - 

T  "  T  V 

^ ref  ^  ref  r  ref 
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P  = 


P 

Pref 


P  = 


P 

Pref 


P  = 


JL 

Pref 


where  the  reference  quantities  are  given  by: 

Lref  is  a  user-defined  constant  length  scale 

pref  is  the  user-defined  total  pressure 

Tref  is  the  user  defined  total  temperature 

Pref  is  the  total  density  defined  as  pref  =  Pre/RrefTref 

1  /7 

aref  is  the  freestream  total  acoustic  velocity  defined  as  aref=  (yRrefTref) 

1/7 

Vref  is  the  reference  velocity  defined  as  are/(y) 

Pref  is  the  reference  viscosity  defined  as  pref  =  PrejVrefLref 

By  rewriting  the  equations  in  dimensionless  form,  it  makes  it  easier  to  compare  with 

experimental  results.  As  a  result,  the  system’s  behavior  and  dependence  on  various 

parameters  can  be  extracted  more  easily. 

Application  of  Harmonic  Balance  Method  to  Governing  Equations 

To  introduce  the  harmonic  balance  method  in  a  simplistic  manner,  one  can 
assume  that  the  flow  is  inviscid  and  non-heat  conducting,  as  was  done  by  Hall,  et  al.  in 
2002  (Hall,  Thomas,  and  Clark,  879-886).  In  his  paper,  the  flow  was  modeled  by  the 
two-dimensional  Euler  equations: 

dU  8F(U)  8G(U )  n 

- +  — - — -  +  — - — -  =  0  ( 

dt  8x  8y 

where  the  vector  of  conservation  variables  and  the  flux  vectors  F  and  G  are  given  by: 
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pv 
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p\>h 

(20) 


To  begin  the  analysis,  one  first  assumes  that  the  flow  behind  the  cylinder  is  temporally 
periodic  whereby, 


U(x,y,t)  =  U(x,y,t  +  T) 


(21) 
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where  T  is  the  period  of  the  unsteadiness  and  is  given  by  T  =  2tt:/co.  Here,  co  is  the 
frequency  of  the  unsteadiness.  As  result  of  the  temporal  periodicity  of  the  flow,  the  flow 
variables  can  be  represented  as  a  Fourier  series  in  time  with  spatially  varying  coefficients. 
Therefore,  the  density  becomes: 

p{x,y,t)  =  Y,RA^yyj(mt  (22) 


Similar  relationships  can  be  obtained  for  the  other  conservation  variables  as  well  as  1/p, 
which  is  necessary  for  determining  the  enthalpy  and  pressure.  In  theory,  one  would  sum 
the  variables  over  all  n  but  in  practice,  the  solution  is  truncated  and  summed  over  a  finite 
number  of  terms.  Next,  the  Fourier  expansions  are  substituted  into  the  Euler  equations 
and  expressions  for  the  conservation  of  mass,  momentum,  and  energy  can  be  developed. 
The  terms  in  the  resulting  equations  are  grouped  by  frequency  and  each  frequency 
component  must  satisfy  the  conservation  equations  individually.  The  equations  can  then 
be  written  in  vector  form  as: 


dImpMi+m= o 

dx  dy 


(23) 


where  the  tilde  indicates  a  Fourier  transformed  variable  and  S  is  the  frequency 
component  resulting  from  the  transfonned  time  dependent  variable  dU/dt.  In  addition, 
since  the  conservation  variables  are  all  real-valued,  it  is  only  necessary  to  store  the 
Fourier  coefficients  for  non-negative  n.  For  example,  if  NH  harmonics  are  retained  in  the 
Fourier  series,  then  2NH  +  1  coefficients  must  be  stored  for  each  flow  variable.  There  is 
one  for  the  zeroth  hannonic  (mean  flow)  and  2NH  for  the  real  and  imaginary  components 
of  the  rest  of  the  harmonics  (Hall,  Thomas,  and  Clark  2002,  881).  However,  this  method 
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proved  to  be  extremely  difficult  and  computationally  expensive  so  an  alternate  approach 
was  developed  to  reduce  the  complexity  of  the  problem. 

Utilization  of  Time  Domain  Methodology 

The  new  method  works  in  terms  of  time  domain  variables  by  assuming  that  the 
values  of  U,  F,  and  G  could  be  evaluated  at  2NH  +  1  equally  spaced  points  in  time  over 
one  temporal  period.  Therefore, 


U*  =  E~lU  (24) 

where  U*  is  a  vector  of  the  conservation  variable  evaluated  at  2Nh  +1  points  and  E1  is 
the  inverse  discrete  Fourier  transfonn  operator.  Similar  expressions  can  be  developed  for 
the  flux  vectors  and  the  new  governing  vector  equation  becomes: 


8F  * 
dx 


dy 


+  S*  =  0 


(25) 


where 


I  ^ 

=  ja>E~lNEU*  * -  (26) 

8t 

S*  is  a  spectral  operator  which  approximates  the  change  in  U*  with  respect  to  time.  This 
method  can  easily  be  extended  to  the  two-dimensional,  unsteady,  viscous  Navier-Stokes 
equations  used  in  this  study.  The  technique  used  to  actually  solve  the  governing 
equations  will  be  described  in  the  next  section. 

Numerical  Solution  Technique 

To  obtain  the  solution  to  the  harmonic  balance  equations,  a  pseudotime  term  was 
introduced.  This  allows  the  equations  to  be  marched  to  a  steady  state  using  conventional 
CFD  techniques.  For  the  Euler  equations,  the  equation  becomes: 
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dU  *  dF  *  dG  * 

- 1 - 1 - 

dr  dx  dy 


+  S*  =  0 


(27) 


where  x  is  a  new  time  variable  merely  used  to  march  the  equation  to  a  steady  state.  This 
equation  very  closely  resembles  the  differential  form  of  the  Euler  equations  so  it  is 
analogous  to  solving  the  time  domain  equations  using  a  CFD  scheme.  Therefore,  the 
overall  method  would  consist  of  pseudotime  marching  N  x  (2Nh  +1)  dependent  variables, 
where  N  is  the  number  of  mesh  points  and  Nh  is  the  number  of  harmonics  retained  in  the 
solution. 


Grid  Generation 


In  general,  the  grid  is  constructed  based  on  the  geometric  attributes  of  the  body 
under  consideration.  The  flow  properties  at  each  point  in  the  grid  are  assigned  a  value 
based  on  the  initial  flow  condition  (Ni  1982,  1565).  Then,  the  governing  equations  are 
discretized  into  distinct  cells  on  a  computational  mesh.  The  solution  is  then  numerically 
integrated  for  each  cell  over  the  entire  domain.  Finally,  the  governing  equations  are 
marched  in  the  psuedo  time  domain.  Multi-grid  and  local  time  stepping  methods  can  be 
used  to  accelerate  convergence.  Effectively,  multi-grid  uses  a  combination  of  fine  grids 
and  coarse  grids  to  achieve  both  rapid  solution  convergence  as  well  as  accuracy. 
Lax-Wendroff  Method 


In  particular,  for  the  two-dimension  Navier-Stokes  equations,  a  computational 
grid  is  generated  for  each  time  level  and  then  the  hannonic  balance  equations  are 
discretized  over  the  entire  domain  using  common  CFD  techniques.  The  conservation 
variables  are  stored  at  each  node  of  the  grid  for  each  time  level.  In  this  analysis,  the  two- 
step  Lax-Wendroff  method  was  used  to  solve  the  harmonic  balance  equations.  The  Lax- 
Wendroff  method  is  a  node-centered  conservative  finite  difference  scheme.  For  the 
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Navier-Stokes  equations,  the  Lax-Wendroff  scheme  is  second-order  accurate  in  space  and 
first  order  accurate  in  time  due  to  viscous  terms  (Davis,  Ni,  and  Carter  1987,  407).  As 
before,  the  Navier-Stokes  equations  are  given  as: 

dU  dF  dG  „ 

- +  —  +  —  =  5  (28) 

dt  ox  8y 

where  U,  F,  G,  and  S  are  defined  in  the  preceding  section.  The  Lax-Wendroff  method  is 
developed  from  a  Taylor  series  expansion  of  U  in  time: 


U  (x,  y,  t  +  At)  =  U  (x,  y )  +  At 


dU(x,y)  |  (  At 2  I  d2U(x,y)  |  ^ 

'1  t  / 


The  first  derivative  tenn  can  be  detennined  from  the  governing  equation  and  therefore 
can  be  written  as: 


dU__  SF+^G  +s 
dt  _dx  dy  _ 


The  second  time  derivative  can  be  found  by  differentiating  this  equation  with  respect  to 


d~U  \  d  (  dFYdU)  d  ( dG\( dU\\  dS\(dU 


dxydU A  dt  J  dy\dU A  dt  )  \dU A  dt 


When  these  expressions  are  substituted  back  into  the  Taylor-series  expansion  for  U,  the 
resulting  equation  becomes: 


TT,  k  \  TT,  \  A(dF  dG\  1  f At2Yo'AF  dAG  A1  3n 

U(x,y,t  +  At)  =  U(x,y)-At  —  +  — - S - —  — —  +  — - AS  +0(Ati)  (32) 

l  ox  uy  J  [21  ux  uy 


where  AU=  (cU/dt)At,  AF  =  (cF/dU)AU,  AG  =  (FG/dU)AU,  AH  =  (dH/oU)AU, 


and  AS  =  (cS/dU)AU. 
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Therefore,  the  first  derivative  representation  provides  a  first-order  correction 
whereas  the  second  derivative  gives  a  second-order  correction.  This  technique  is  applied 
to  the  hannonic  balance  described  in  the  previous  section  and  the  solution  can  be 
obtained  by  marching  these  equations  in  time  at  each  grid  point.  A  steady  state  is 
reached  when  the  corrections  are  driven  to  zero.  Since  this  method  was  originally 
developed  for  inviscid  flow  analysis,  some  small  changes  must  be  made  to  the  integration 
scheme  to  handle  the  viscous  components  (Davis,  Ni,  and  Carter,  409).  Finally,  it  is 
common  to  use  an  artificial  viscosity  to  maintain  stability  in  regions  with  large 
discontinuities  and  coarse  grid  spacing.  In  this  study,  both  second  and  fourth  difference 
numerical  smoothing  are  used  to  eliminate  oscillations. 

Application  of  Initial  and  Boundary  Conditions 

In  addition,  proper  boundary  and  initial  conditions  must  be  applied.  In  particular, 
at  each  iteration  of  the  code,  the  no-slip  boundary  condition  must  be  applied  on  the 
cylinder  surface  since  viscous  flow  analysis  is  being  conducted.  Furthennore,  the 
periodicity  conditions  must  be  implemented.  Also,  non-reflecting  far  field  boundary 
conditions  should  be  applied  in  the  frequency  domain  such  that  there  are  no  disturbances 
far  from  the  body  (Hall,  Thomas,  and  Clark  2002,  882).  Finally,  both  upstream  and 
downstream  initial  flow  conditions  are  prescribed. 

Frequency  Search  Procedure 

Once  the  system  has  been  completely  described,  CFD  techniques  can  be  invoked 
and  the  equations  are  marched  in  pseudo-time  until  convergence  is  reached.  However, 
nonlinear  harmonic  balance  analysis  requires  the  user  to  input  the  frequency  of  the  flow 
instability  (a  known  quantity)  to  the  system.  Since  the  purpose  of  the  study  is  to  find  the 
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NSV  frequency,  two  different  approaches  can  be  applied.  A  simple  method  is  to  choose 
a  wide  range  of  frequencies  and  compute  the  solution  residual  of  the  hannonic  balance 
equation  for  each  case.  The  NSV  frequency  corresponds  to  that  for  which  the  harmonic 
balance  solution  residual  goes  to  zero.  However,  previous  studies  by  McMullen,  et  al. 
indicate  the  solution  residual  drops  off  suddenly  at  the  precise  NSV  frequency 
(McMullen,  Jameson,  and  Alonso  2002,  6).  As  a  result,  this  method  of  guessing 
frequencies  can  easily  miss  the  correct  NSV  frequency  and  also  requires  numerous 
iterations  of  the  harmonic  balance  code. 

Phase  Error  Method 

Therefore,  a  novel  frequency  search  technique  was  developed  by  Hall  (private 
communication).  This  technique  requires  much  fewer  hannonic  balance  calculations 
than  the  frequency  sweep  method.  For  a  given  frequency,  the  converged  solution  will 
have  a  constant  amplitude  and  phase.  The  phase  shift  can  be  calculated  relatively  simply 
once  the  unsteady  lift  is  detennined.  After  each  iteration  of  the  HB  solver,  the  phase  is 
updated  and  the  change  in  phase  from  one  iteration  to  the  next  is  computed.  When  an 
incorrect  frequency  is  inputted  using  the  nonlinear  hannonic  balance  method,  it  was 
discovered  that  the  amplitude  of  any  unsteady  first  harmonic  quantity,  such  as  unsteady 
lift,  stays  nearly  constant  but  the  phases  rotates  from  one  iteration  to  the  next  and 
eventually  settles  down  to  constant  value.  Once  this  value  is  calculated  for  a  couple 
different  frequencies,  it  becomes  readily  apparent  that  the  phase  is  nearly  linearly  related 
to  the  frequency.  The  exact  frequency  can  be  calculated  by  determining  the  frequency 
for  which  the  phase  shift  is  zero.  This  can  be  accomplished  by  simple  interpolation. 
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Thus,  this  simple  phase  error  technique  results  in  a  reasonable  approximation  of  the  NSV 
frequency  and  only  requires  a  few  harmonic  balance  calculations. 

3.  APPLICATIONS  OF  HARMONIC  BALANCE  METHOD 

3.1  STATIONARY  CYLINDER  IN  CROSS  FLOW 

A  well-studied  fluid  dynamics  phenomenon  is  a  stationary  cylindrical  structure  in 
cross  flow.  This  test  case  is  important  for  the  study  of  towers,  cables,  antennae,  wires, 
etc.  As  a  result,  extensive  experimental  data  is  available  demonstrating  the  relationship 
between  Strouhal  number  and  Reynolds  number.  Therefore,  as  an  initial  illustration  of 
the  merits  of  the  hannonic  balance  method,  the  first  goal  was  to  reproduce  these 
experimentally  determined  results.  The  problem  was  modeled  using  the  HB  method  and 
a  129x65  mesh  for  various  Reynolds  numbers.  A  typical  mesh  is  shown  in  Figure  3  with 
129  points  in  the  circumferential  direction  and  65  points  in  the  radial  direction.  The  mesh 
boundary  is  40  diameters  from  the  center  of  the  cylinder  in  this  figure. 


Figure  3.  Computational  Grid  (129x65  mesh)  for  Cylinder  in  Cross  Flow 
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Application  of  Phase  Error  Method 

To  determine  the  Strouhal  number  representing  the  NSV  frequency  for  a 
particular  Reynolds  number,  it  was  necessary  to  run  the  code  for  different  Strouhal 
numbers  and  note  the  solution  residual  for  each  case.  Figure  4  shows  the  HB  solution 
residuals  for  a  number  of  assumed  Strouhal  numbers  at  a  Reynolds  number  of  170  using 
two  harmonics.  As  can  be  seen  from  the  graph,  the  predicted  Strouhal  number  is 
approximately  0.1865.  The  NSV  frequency  or  the  “correct”  frequency  is  that  for  which 
the  solution  residual  drops  down  significantly  i.e.  |dq|  =  |0|  (the  difference  between  the 
current  solution  and  the  previous  solution  is  zero).  Flowever,  this  method  requires  the 
user  to  run  the  code  for  many  different  cases  to  hone  in  on  the  precise  NSV  frequency. 


Figure  4.  HB  Solution  Residual  versus  Strouhal  Number  -  Re  =  170  (2  Harmonics) 

A  more  efficient  method  is  to  use  the  phase  error  method  described  previously 
that  utilizes  the  change  in  an  integrated  global  quantity  such  as  the  phase  of  the  unsteady 
lift  as  opposed  to  the  solution  residual.  By  choosing  the  same  Strouhal  numbers,  the 
phase  shift  per  iteration  was  also  determined.  A  plot  of  these  quantities  is  shown  in 
Figure  5  and  it  can  be  seen  that  the  change  in  phase  of  the  unsteady  lift  is  nearly  linear 
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with  respect  to  the  Strouhal  number.  The  zero  crossing  represents  the  NSV  non- 
dimensional  frequency  and  is  also  the  point  where  the  FIB  solution  residual  exhibits  a 
sudden  drop-off.  Therefore,  it  is  only  necessary  to  run  the  code  for  a  small  number  of 
different  frequencies  and  note  the  change  in  phase  of  the  unsteady  lift  for  each  case. 

Then,  one  can  interpolate  between  these  points  to  find  the  zero  crossing,  which 
corresponds  to  the  precise  NSV  frequency.  Therefore,  this  frequency  search  technique 
provides  a  solution  much  more  quickly  than  simply  searching  directly  for  the  frequency 
where  the  FIB  solution  residual  is  a  minimum.  Furthermore,  a  more  accurate  solution  can 
be  obtained  by  adding  more  harmonics  and  performing  a  mesh  refinement  study. 


Strouhal  Number,  St  =  fD/U^ 

Figure  5.  Change  in  Phase  Angle  of  Unsteady  Lift  Per  Iteration  -  Re  =  170 

(2  Harmonics) 

Strouhal  -  Reynolds  Number  Relationship 

This  process  was  repeated  for  a  range  of  Reynolds  numbers  in  the  laminar  flow 
regime  and  a  corresponding  Strouhal  number  was  detennined  for  each  case.  The 
resulting  Strouhal  numbers  were  then  compared  with  experimental  data  collected  by 
Williamson  in  1996  (Williamson,  494).  The  results  can  be  found  in  Figure  6. 


Duke  University 


Page  26 


March  2006 


F49620-03- 1-0204 


Final  Report 


Reynolds  Number,  Re 

Figure  6.  Cylinder  in  Cross  Flow  -  Strouhal  Number  Dependency  on  Reynolds 

Number 

As  can  be  seen  from  the  plot,  the  FIB  method  shows  good  agreement  with  the 
experimentally  determined  values.  Previous  studies  had  indicated  large  discrepancies  in 
the  St  -Re  relationship  in  the  laminar  flow  region  but  Williamson  found  that  the  wake 
frequency  was  very  sensitive  to  the  experimental  setup,  particularly  the  3-D  effects  of 
oblique  vortex  shedding.  To  counter  these  effects,  parallel  shedding  was  induced  by 
angling  the  endplates  on  the  cylinder  (Williamson  1988,  2744).  Therefore,  by  applying 
this  end  boundary  condition  and  ensuring  the  freestream  is  sufficiently  uniform,  a  curve 
is  produced  that  is  universal  as  well  as  continuous.  Therefore,  the  harmonic  balance 
method  accurately  predicts  the  relationship  between  Strouhal  number  and  Reynolds 
number. 
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Stationary  Cylinder  Flow  Simulation 

A  typical  flow  solution  for  the  total  pressure  contours  for  Reynolds  numbers  of 
50,  100,  and  170  are  shown  in  Figures  7,  8,  and  9,  respectively.  These  plots  demonstrate 
the  low-pressure  region  that  is  generated  behind  the  cylinder  and  also  show  the  formation 
of  vortices  in  the  cylinder’s  wake.  In  particular,  for  Re  =  50,  it  shows  the  beginning  of 
the  vortices  being  shed  alternatively  from  the  top  and  bottom  of  the  cylinder.  For 
Reynolds  numbers  of  100  and  170,  the  wake  has  become  unstable  and  the  vortices  have 
detached  themselves  from  the  cylinder.  Consequently,  staggered  rows  of  vortices  of 
opposite  sign  are  being  formed.  Therefore,  the  HB  method  is  able  to  produce  time 
accurate  results  and  is  a  viable  tool  to  predict  the  frequency  of  the  unsteady  flow  about  a 
stationary  cylinder  in  cross  flow  and  the  solution  at  chosen  Reynolds  numbers  agrees 
well  with  previous  experimental  investigations.  Furthermore,  the  observed  flow  patterns, 
such  as  the  Von  Karman  vortex  street,  are  consistent  with  previous  experimental  flow 
simulation  results.  This  can  be  seen  by  comparing  the  figures  below  with  Williamson’s 
photograph  of  the  formation  of  Von  Karman  vortex  streets  in  the  introduction. 
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Figure  7.  Total  Pressure  Contours  for  Flow  Over  a  Cylinder  at  Re  =  50 
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Figure  8.  Total  Pressure  Contours  for  Flow  Over  a  Cylinder  at  Re  =  100 
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Figure  9.  Total  Pressure  Contours  for  Flow  Over  a  Cylinder  at  Re  =  170 
Unsteady  Cylinder  Lift 

In  addition  to  determining  the  frequency  of  the  fluid  dynamic  instability,  another 
important  quantity  to  be  able  to  calculate  is  the  RMS  amplitude  of  the  lift  on  the  cylinder. 
At  the  NSV  frequency  for  a  range  of  Reynolds  numbers,  the  amplitude  of  the  first 
harmonic  lift  acting  on  the  shedding  cylinder  is  determined.  A  plot  demonstrating  this 
relationship  can  be  found  in  Figure  10.  This  plot  gives  a  measure  of  the  unsteady  lift  (1st 
harmonic).  The  alternate  periodic  shedding  causes  the  pressure  fluctuations  at  around  fst 
to  be  essentially  out-of-phase  between  the  upper  and  lower  side  of  the  cylinder  so  the  lift 
fluctuation  energy  is  concentrated  to  a  band  around  fst  (Norberg  2003,  57).  Furthermore, 
by  extrapolating  to  the  Reynolds  number  of  zero  oscillating  lift,  the  onset  of  the  vortex 
shedding  can  be  determined.  Figure  10  shows  that  this  occurs  at  about  a  Reynolds 
number  of  47,  which  is  approximately  the  same  as  the  value  determined  from  nonlinear 
dynamic  numerical  techniques  (Norberg  2001,  464).  In  addition,  the  plot  shows  that 
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there  is  rapid  increase  in  the  unsteady  RMS  lift  coefficient  within  the  laminar  shedding 

regime.  This  is  characteristic  of  all  supercritical  Flopf  bifurcations  in  which  the  size  of 

1  /2 

the  limit  cycle  grows  continuously  from  zero,  and  increases  proportional  to  (ju  -  pc) 
where  p  is  a  control  parameter.  In  fact,  experimental  studies  and  2-D  simulations  have 
shown  that  the  unsteady  lift  coefficient  exhibits  this  square-root  dependency  on  Reynolds 
number,  Cl’  oc  e  =  (Re-Rec)/Rec))  where  Rec  is  the  critical  Reynolds  number  for  which 

shedding  begins  (Norberg  2001,  464).  Therefore,  by  examining  the  nonlinear  dynamics 
of  the  system,  the  unsteady  lift  provides  a  measure  of  the  stability  of  the  system  and  can 
be  used  to  determine  the  transition  from  steady  flow  to  Von  Karman  vortex  shedding  in 
the  2-D  cylinder  wake.  Furthermore,  the  HB  method  demonstrates  the  unsteady  lift’s 
square-root  dependency  on  Reynolds  number. 


Reynolds  Number,  Re 

Figure  10.  Computed  Amplitude  of  First  Harmonic  of  Unsteady  Lift  Acting  on 
Shedding  Cylinder  from  Onset  at  Re  =  47 
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Experimental  Procedure  and  Results 

To  validate  the  results  obtained  by  the  HB  method,  a  comparison  was  made 
with  the  experimental  study  of  Tanida,  et  al.  as  well  as  numerous  2-D  numerical 
simulations.  Figure  1 1  shows  that  the  results  obtained  from  the  HB  method  show 
reasonable  agreement  with  the  2-D  computational  data  but  vary  considerably  from  the 
experimental  data.  This  discrepancy  may  be  explained  in  part  by  the  use  of  unsealed 
gaps  in  Tanida’s  towing  tank  study  (Norberg  2003,  65).  Tanida,  et  al.  utilized  a  force 
element  method  to  measure  the  fluctuating  lift  in  which  the  load-transmitting  part  of  the 
cylinder  is  connected  to  a  cantilever  beam  element  that  is  fixed  to  a  base  inside  a 
“dummy”  part  of  the  cylinder.  Keefe  found  that  unsealed  gaps  can  result  in  a  drastic 
reduction  in  fluctuating  unsteady  lift  forces  as  Reynolds  number  is  decreased,  as  much  as 
10  times  lower  than  with  sealing  (Keefe  1962,  1712).  However,  his  experiment  was 
conducted  at  higher  Reynolds  number  so  it  is  difficult  to  detennine  a  definitive  reason  for 
the  incongruity  between  the  experimental  and  computational  data.  Unfortunately,  this  is 
the  only  experimental  data  available  in  the  Reynolds  regime  of  interest  so  it  is  hard  to 
make  a  qualitative  comparison.  When  the  HB  method  is  compared  with  other  numerical 
simulations,  it  seems  to  slightly  over  predict  the  unsteady  lift  coefficient.  As  a  result, 
further  study  is  required  to  detennine  the  validity  of  Tanida’s  experimental  results. 
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Figure  11.  Comparison  of  Lift  Data  with  Other  Numerical  Results  and 

Experimental  Data 

Mean  Cylinder  Drag 

Another  interesting  feature  to  examine  is  the  response  of  the  cylinder  drag  to  the 
onset  of  unsteadiness  and  vortex  shedding.  There  are  numerous  contributions  to  the  total 
drag  on  a  cylinder,  including  the  viscous  drag  coefficient  and  the  pressure  drag 
coefficient.  Henderson  computed  these  coefficients  numerically  using  a  highly  accurate 
spectral  element  method  based  on  8th  order  polynomials  (Henderson  1995,  2102). 
Furthermore,  McMullen  also  computed  the  drag  forces  as  a  function  of  Reynolds  number 
using  a  frequency  domain  technique.  The  change  from  a  steady  to  an  unsteady  wake  is 
marked  by  a  gradual  decrease  in  the  viscous  drag  coefficient  as  well  as  the  total  drag 
throughout  the  laminar  flow  regime  (Henderson  1995,  2103).  A  plot  comparing  the 
values  obtained  by  Henderson,  McMullen,  and  the  harmonic  balance  method  is  found  in 
Figure  12.  The  results  from  the  harmonic  balance  method  are  very  similar  to  those 
obtained  by  both  Henderson  and  McMullen.  It  is  noted  that  the  mean  drag  coefficients 
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are  relatively  constant  with  respect  to  Reynolds  number  as  compared  to  the  lift 
coefficients.  Therefore,  the  harmonic  balance  method  is  further  validated  as  a  valuable 
tool  for  predicting  the  forces  on  a  cylinder. 
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Figure  12,  Comparison  of  Mean  Coefficient  of  Drag  versus  Reynolds  Number  Data 
with  Henderson  and  McMullen  for  Laminar  Vortex  Shedding 

3.2  CYLINDER  WITH  PRESCRIBED  MOTION 


Next,  the  cylinder  is  forced  to  oscillate  at  a  specified  amplitude  and  frequency.  In 
this  case,  there  is  a  range  of  frequencies  over  which  the  shedding  frequency  will  “lock- 
in”  to  the  frequency  of  the  vibrating  cylinder.  This  synchronization  effect  was  first 
observed  by  Bishop  and  Hassan  and  later  measured  by  Koopman  at  low  Reynolds 
numbers  (Koopman  1967,  508).  Cylinder  vibration  with  frequencies  near  the  shedding 
frequency  can  influence  both  the  pattern  and  the  phasing  of  vortices.  Outside  of  this 
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region,  the  cylinder  will  oscillate  at  a  frequency  close  to  that  of  the  stationary  shedding 
frequency.  Forced  oscillations  were  achieved  experimentally  by  using  a  shaking 
mechanism  to  generate  a  controlled  motion  at  a  range  of  amplitudes  and  Strouhal 
frequencies.  Therefore,  as  the  driving  frequency  is  increased,  two  distinct  behaviors  are 
observed,  inside  and  outside  of  the  so-called  lock-in  regime. 

Technique  Used  to  Determine  Bounds  of  Lock-In  Region 

To  investigate  this  phenomenon  numerically,  the  same  FIB  solution  methodology 
was  used  but  with  the  cylinder  vibrating  in  the  transverse  direction  at  a  prescribed 
frequency.  Koopman  experimentally  determined  the  lock-in  region  for  Reynolds 
numbers  of  100  and  200.  In  a  first  attempt  to  replicate  the  data,  the  code  was  run  at  a 
Reynolds  number  of  150  with  a  fixed  amplitude  for  many  different  Strouhal  numbers, 
and  the  solution  residual  was  noted  for  each  case.  The  solution  residual  within  the  lock- 
in  region  will  converge  to  machine  zero.  Outside  of  this  region,  two  distinct  frequencies 
are  observed  experimentally;  however  the  current  FIB  method  can  only  handle  one 
frequency.  Thus  in  the  region  where  two  frequencies  are  expected  to  be  present,  the 
current  numerical  method  fails  to  converge.  Extension  of  the  HB  method  to  treat  two  or 
more  frequencies  is  possible  and  is  the  subject  of  current  research.  Figure  13  shows  this 
behavior  for  two  different  Strouhal  numbers,  one  inside  and  one  outside  of  the 
synchronization  region. 
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Iteration,  N 

Figure  13.  Behavior  of  Solution  Residual  In  and  Out  of  the  Lock-in  Region 

Therefore,  by  repeating  this  procedure  for  various  amplitudes,  it  was  possible  to  establish 
an  estimate  of  the  left  and  right  bounds  on  the  lock-in  region  based  on  the  FIB  solution 
behavior.  Flowever,  it  is  difficult  to  determine  the  exact  left  and  right  bounds  of  the  lock- 
in  region  due  to  the  inability  of  the  HB  method  to  accurately  determine  the  frequencies 
outside  of  the  lock-in  region  where  the  cylinder  is  shedding  at  a  different  frequency  than 
it  is  vibrating.  Therefore,  there  is  some  subjectivity  in  assessing  what  constitutes  a  fully 
converged  solution  within  the  lock-in  region. 

Initially,  a  mesh  size  of  129x65  was  used  with  a  mesh  boundary  radius  of  forty 
diameters  measured  from  the  cylinder’s  center.  A  comparison  between  Blevins’  data  and 
the  results  obtained  by  the  FIB  method  is  shown  in  Figure  14.  The  dimensionless 
amplitude  is  plotted  versus  the  ratio  of  the  frequency  of  the  vibrating  cylinder  to  the 
stationary  shedding  frequency.  From  the  plot,  it  can  be  seen  that  there  is  good  agreement 
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up  to  an  oscillation  amplitude  of  0.05.  However,  above  this  amplitude,  the  HB  solution 
begins  to  deviate  greatly  from  the  experimental  values.  This  discrepancy  may  be  due  to 
the  use  of  a  coarse  mesh,  a  large  mesh  boundary,  or  the  code’s  limitation  of  the  use  of 
only  3  harmonics.  As  a  result,  a  number  of  changes  were  made  to  achieve  better 
agreement  with  Blevins’  experimental  data. 
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- Re  =  100  Experimental  Re  =  200  Experimental 

— A —  Re  =  150  (129x65  mesh,  NH  =3)  •  Re  =  100  (193x97  mesh,  NH  =  3) 

X  Re  =  100  (193x97  mesh,  NH  =  7) 

Figure  14.  Demonstration  of  Lock-in  Phenomenon 
Effect  of  Mesh  Refinement 

The  initial  improvements  involved  utilization  of  finer  meshes  and  a  smaller  mesh 
boundary.  Both  a  193x97  and  a  257x129  mesh  were  considered  as  well  as  a  mesh 
boundary  of  20  times  the  diameter  of  the  cylinder  (20D)  instead  of  40D.  Table  1  shows 
the  results  obtained  from  the  mesh  refinement  study  for  dimensionless  amplitude  of  0. 10. 
As  can  be  seen,  the  use  of  a  finer  mesh  improves  the  results  at  an  amplitude  of  0. 10  but 
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there  is  still  a  significant  difference  in  the  bounds.  Therefore,  a  further  refinement  is 
required  to  achieve  better  accuracy  at  higher  vibratory  amplitudes. 

Addition  of  Filter  to  Reduce  Far  Field  Boundary  Effects 

Due  to  the  constraints  of  using  the  HB  method  with  a  large  grid  domain,  the 
previous  code  limited  the  user  to  only  three  harmonics.  In  an  attempt  to  avoid  this 
problem,  a  subroutine  was  added  to  act  as  a  filter  for  the  HB  method.  The  filter  was 
designed  to  zero  out  the  higher  harmonics  as  you  go  farther  and  farther  out  in  the  domain, 
away  from  the  cylinder  where  the  wake  effects  are  small.  Therefore,  it  is  possible  to  use 
as  many  as  seven  harmonics  to  achieve  a  more  accurate  solution.  By  examining  Table  1, 
the  results  show  that  adding  up  to  seven  harmonics  greatly  improves  the  estimate  of  the 
left  bound  of  the  lock-in  region  and  showed  a  slight  improvement  in  the  right  bound.  In 
addition,  the  results  from  the  193x97  mesh  were  added  to  Figure  20  to  demonstrate  the 
new  bounds.  Therefore,  the  filter  was  a  valuable  tool  because  it  gave  a  way  to  study  the 
effect  of  adding  more  harmonics  to  the  solution.  It  was  found  that  the  use  of  more 
harmonics  resulted  in  a  better  approximation  to  the  experimental  values  obtained  by 
Blevins.  However,  mesh  size  had  only  a  small  effect  on  the  bounds  of  the  lock-in  region. 


Left  Bound 

Right  Bound 

Mesh  Size 

Re 

h/D 

Radius 

NH 

f/fN 

f/fN 

129x65 

100 

0.10 

20 

3 

0.7656 

1.1016 

193x97 

100 

0.10 

20 

3 

0.7656 

1.1016 

257x129 

100 

0.10 

20 

3 

0.7891 

1.1094 

193x97 

100 

0.10 

20 

7 

0.8438 

1.1094 

Experimental 

100 

0.10 

- 

- 

0.875 

1.1284 

129x65 

150 

0.10 

40 

3 

0.7708 

1 .0972 

Table  1.  Effect  of  Mesh  and  Harmonic  Refinement  on  Modeling  Blevin’s 
Experimental  Lock-in  Region  for  Re  =  100 
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Unsteady  Lift  for  Prescribed  Motion 

In  addition  to  determining  the  lock-in  region,  the  unsteady  lift  was  calculated. 
Tanida,  et  al.  experimentally  measured  this  quantity  in  a  towing  tank  with  30  mm 
diameter  test  cylinders  in  which  the  lift  and  drag  forces  are  sensed  by  strain  gauges 
(Tanida  1973,  773).  Oil  was  used  as  the  fluid  because  it  allows  the  unsteady  aerodynamic 
forces  on  the  oscillating  cylinder  to  be  measured  with  reasonable  accuracy  (Tanida  1973, 
771).  The  study  was  conducted  at  a  Reynolds  number  of  80  and  a  non-dimensional 
amplitude  of  0.14.  Once  again,  the  HB  method  was  used  at  these  conditions  and 
compared  with  the  experimental  data.  A  plot  of  the  results  can  be  found  in  Figure  16. 

The  HB  method  shows  remarkable  agreement  with  the  experimental  results  of  Tanida,  et 
al.  The  stability  of  the  cylinder  oscillation  is  dependent  on  the  component  of  the 
unsteady  lift  force  that  is  in  phase  with  the  oscillating  velocity  (Tanida  1973,  774). 
Therefore,  the  cylinder  will  be  unstable  if  the  imaginary  part  of  the  unsteady  lift 
coefficient  is  greater  than  zero,  which  is  characteristic  of  negative  aerodynamic  damping. 
Real  and  Imaginary  Part  of  the  Unsteady  Lift 

The  real  and  imaginary  components  of  the  lift  coefficient  are  plotted  below.  The 
plot  of  the  imaginary  component  shows  the  system  is  stable  for  St  =  0.1000  to  St  = 

0.1300  and  then  becomes  unstable  throughout  the  rest  of  the  lock-in  region.  The  stability 
of  the  system  is  determined  by  assuming  harmonic  motion  in  which 

CL  =|  CL  |  e-ia *  =|  CL  |  (33) 

where  a  negative  value  indicates  a  stable  solution  and  a  positive  value  gives  an 
exponentially  growing  term  which  causes  instability.  A  sensitivity  study  was  conducted 
to  determine  the  effect  of  including  more  and  more  harmonics.  There  does  not  appear  to 
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be  a  substantial  benefit  to  keeping  more  than  two  harmonics  in  the  HB  solution  for  the 
parameter  range  studied  here. 


Figure  15.  Magnitude  of  the  Real  Part  of  the  Unsteady  Lift  Coefficient  versus 
Strouhal  Number  for  a  Single  Cylinder  Oscillating  Transversely  (h/D  =  0.14) 


Figure  16.  Magnitude  of  the  Imaginary  Part  of  the  Unsteady  Lift  Coefficient  versus 
Strouhal  Number  for  a  Single  Cylinder  Oscillating  Transversely  (h/D  =  0.14) 
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Amplitude  of  the  Unsteady  Lift 

In  addition,  a  plot  of  the  amplitude  of  the  unsteady  lift  coefficient  was 
constructed.  By  examining  the  plot  in  Figure  17,  the  amplitude  appears  to  steadily 
increase  throughout  the  lock-in  region  until  a  peak  displacement  is  reached.  When 
compared  to  the  stationary  cylinder  case,  nearly  the  same  Cl  value  is  obtained  for  Re  = 
80.  In  addition,  the  cylinder  was  also  forced  to  oscillate  at  half  the  vibratory  amplitude, 
i.e.  h/D  =  0.07  and  it  was  discovered  that  lock-in  region  became  smaller  but  the 
amplitude  of  the  lift  coefficient  fell  on  the  same  curve  as  the  h=0.14  case.  The  cylinder 
was  also  forced  to  vibrate  at  higher  amplitudes  and  the  lift  coefficient  actually  decreased 
(see  Figure  18).  Therefore,  cylinder  motion  does  not  dramatically  affect  the  lift 
coefficient  for  oscillation  amplitudes  up  to  about  h/D  =  0.40.  Therefore,  the  lift 
coefficient  appears  to  be  relatively  independent  of  the  prescribed  amplitude  when  the 
cylinder  is  driven  at  an  amplitude  of  approximately  10  percent  or  less  of  the  cylinder’s 
diameter  and  for  higher  vibratory  amplitudes,  the  lift  coefficient  decreases.  This  has 
important  implications  for  the  study  of  non-synchronous  vibrations  of  turbomachinery 
blades  because  it  may  not  be  necessary  to  couple  the  NSV  fluid  dynamic  solution  with 
blade  motion  for  sufficiently  small  blade  amplitudes,  which  is  a  much  easier 
computation. 
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■4—  h/D  =  0.10,  NH  =  2,  Re  =  150  h/D  =  0,  NH  =  2  (Stationary  Cylinder) 

Figure  17.  Amplitude  of  the  Unsteady  Lift  versus  Strouhal  Number  for  a  Reynolds 
Number  of  80  and  a  Mesh  Size  of  129x65  with  a  Radius  of  20 
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Figure  18.  Effect  of  Oscillation  Amplitude  on  the  Magnitude  of  the  Unsteady  Lift 
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Phase  Shift  as  a  Function  of  Strouhal  Number 

Another  important  characteristic  of  the  lock-in  region  is  the  phase  angle  as  a 
function  of  the  Strouhal  number  within  the  lock-in  region.  The  phase  angle  is  a  measure 
of  the  phase  difference  between  the  displacement  of  the  cylinder  and  the  lift  force.  By 
examining  Figure  18  below,  the  cylinder  lift  undergoes  an  abrupt  180  degrees  phase  shift 
between  the  shedding  and  the  cylinder  motion  at  about  a  Strouhal  number  of  0. 1500.  For 
a  Reynolds  number  of  about  80,  the  stationary  cylinder  shedding  frequency  was  found  to 
be  about  St  =  0.1450  using  the  HB  method.  Therefore,  it  appears  that  this  180  degrees 
phase  shift  occurs  as  the  cylinder  vibration  frequency  passes  through  the  natural  shedding 
frequency.  This  could  be  explained  by  Zdravkovich’s  physical  observation  of  the  flow 
behind  the  cylinder.  By  examining  the  flow  over  an  oscillating  cylinder,  it  is  noted  that 
the  vortices  tend  to  shed  when  the  cylinder  is  near  its  maximum  displacement  (Blevins 
1990,  56).  Zdravkovich  found  that  for  frequencies  below  the  natural  shedding  frequency, 
the  vortex  is  shed  from  the  side  opposite  that  experiencing  its  maximum  displacement. 
However,  for  frequencies  above  the  shedding  frequency,  the  vortex  is  shed  from  the  same 
side  as  the  max  displacement  (Sarpkaya  2003,  64).  Therefore,  the  phase  shift  may  be  due 
in  part  to  a  switch  in  the  timing  of  the  shedding  of  the  vortices.  The  stability  of  the 
system  is  consistent  with  the  plot  of  the  imaginary  coefficient  of  lift.  A  phase  shift  of 
zero  degrees  refers  to  the  condition  in  which  the  force  and  displacement  are  in-phase  and 
no  work  is  being  done. 
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Figure  19.  Phase  Angle  Between  the  Lift  and  Displacement  Within  the  Lock-In 

Region 


Aerodynamic  Damping 

In  addition,  the  aerodynamic  damping  was  also  calculated.  Aerodynamic 
damping  is  a  result  of  fluid  forces  acting  on  the  structure.  As  can  be  seen  from  Figure  19, 
the  damping  becomes  increasingly  negative  until  it  reaches  a  minimum  value  and  then 
starts  to  increase  again.  A  negative  aerodynamic  damping  results  in  a  net  energy  input  to 
the  cylinder  vibration.  The  most  unstable  damping  coefficient  occurs  at  a  Strouhal 
number  of  0.1600  which  coincides  with  the  point  for  which  maximum  work  is  being 
done. 
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Figure  20.  Aerodynamic  Damping  in  the  Lock-in  Region  as  a  Function  of  Strouhal 

Number 

Therefore,  by  examining  the  amplitude,  phase  angle,  and  damping  coefficient  within  the 
lock-in  region,  various  flow  characteristics  were  deduced.  Specifically,  the  stability  of 
the  flow,  the  frequency  for  which  the  maximum  response  occurs,  and  the  behavior  at  the 
endpoints  of  the  lock-in  region  were  found  and  compared  with  the  physics  of  the  flow 
behind  a  cylinder. 

Mean  Cylinder  Drag  for  Prescribed  Motion 

In  addition  to  the  fluctuating  lift  forces  on  the  cylinder,  the  mean  drag  was  also 
calculated  as  a  function  of  Strouhal  number.  Figure  20  shows  the  results  obtained  by  the 
HB  method  as  compared  to  the  experimental  results  of  Tanida.  By  examining  the  plot,  it 
can  be  seen  that  the  drag  coefficient  steadily  increases  within  the  lock-in  region  until  it 
reaches  a  maximum  around  Cdo  =  1 .63  at  St  =  0. 1650  for  the  FIB  method.  On  the  other 
hand,  Tanida ’s  results  indicate  a  maximum  value  of  Cdo  =  1.87  at  St  =  0.1500.  Flowever, 
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despite  this  discrepancy,  the  overall  trend  appears  to  be  the  same  and  the  HB  method 
shows  relatively  good  agreement  with  Tanida’s  experimental  data.  Therefore,  the  FIB 
method  is  also  capable  of  determining  the  drag  forces  on  a  cylinder  when  it  is  subjected 
to  a  prescribed  motion. 


♦  Total  Drag  -  NH  =  2 
-■ —  Total  Drag  -  NH  =  5 
Total  Drag  -  NH  =  7 
Viscous  Drag  -  NH  =  2 
—  Viscous  Drag  -  NH  =  5 
-• —  Viscous  Drag  -  NH  =  7 
—I —  Experiment  -  Tanida  (1973) 


0.0500  0.1000  0.1500 

Strouhal  Number,  St 


0.2000 


Figure  21.  Comparison  of  the  Mean  Coefficient  of  Drag  versus  Strouhal  Number 

Within  the  Lock-in  Region 


3.3  ELASTICALLY  MOUNTED  CYLINDER 


Finally,  it  is  important  to  determine  the  response  characteristics  of  a  vortex- 
excited  spring  supported  cylinder  in  the  stable,  laminar  flow  regime.  In  this  case,  as  the 
flow  velocity  is  increased  or  decreased,  the  shedding  frequency  can  approach  the  natural 
frequency  of  the  structure.  At  a  critical  velocity,  the  shedding  frequency  will  lock-in  to 
the  structure  frequency  (Blevins  1973,  21).  In  the  synchronization  region,  resonant 
oscillation  conditions  can  occur  and  produce  large  amplitude  responses.  Therefore,  by 
determining  the  interaction  between  the  flow  oscillations  and  the  cylinder  motion,  it  will 
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be  possible  to  predict  the  cylinder  oscillation  frequency  and  the  resultant  cylinder 
amplitude. 

Aeroelastic  Cylinder  Model 

This  cylinder  aeroelastic  “lock-in”  problem  is  modeled  by  a  single  degree  of 

freedom  spring -mass-damper  system  excited  by  a  transverse  lift  force  expressed  as  L',  the 

lift  per  unit  span.  A  schematic  of  this  system  can  be  seen  in  Figure  22.  As  can  be  seen, 

the  cylinder  is  mounted  on  a  linear  spring  of  stiffness  of  k  and  a  structural  damping 

coefficient  of  d.  The  external  force  on  the  cylinder  is  represented  by  the  unsteady  lift 

generated  by  the  trailing  vortices.  The  governing  equation  that  describes  this  system  is: 

d2h  ,  dh  7 ,  1  TT  ^  _  .... 

m  —  +  d  —  +  kh  =  -  pJ3mDsCv  (34) 

dt  dt  2 

where  m  is  the  mass  of  the  cylinder,  d  is  the  damping  coefficient,  k  is  the  cylinder  spring 
stiffness,  px is  the  fluid  density,  Uao  is  the  fluid  velocity,  D  is  the  cylinder  diameter,  5  is 
the  cylinder  span,  Cl> is  the  cylinder  lift  coefficient,  and  h  is  the  transverse  displacement 
of  the  cylinder. 


Figure  22.  One-Degree-of-Freedom  Model  of  a  Vortex-Excited  Cylinder 
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Next,  one  can  assume  simple  harmonic  motion  so  h  and  Cl  •  can  be  represented  as 

h  =  h,ej0\  CL,  =CLeJ0,t  (35) 

As  a  result  of  this  substitution,  equation  (38)  becomes: 

9  1  2 

{-co  in  +  jcod  +  k)h\  ~—pJJxDsCL.  =  0  (36) 

Non-Dimensional  Parameters 

In  practice,  it  is  more  common  to  study  the  non-dimensionalized  form  of  this 
equation  because  it  reduces  the  number  of  parameters  and  also,  makes  it  easier  to 
compare  with  experimental  results.  As  a  result,  some  new  dimensionless  parameters 
should  be  introduced.  The  first  parameter  is  the  natural  frequency  of  the  system,  which  is 
independent  of  any  initial  excitations.  It  is  expressed  as 

(37) 

Another  parameter  of  interest  related  to  the  natural  frequency  is  the  damping  factor.  This 
value  characterizes  the  energy  dissipated  by  the  cylinder  as  it  vibrates.  It  is  described  by 
the  following  equation: 

<r=A-  os) 

2  mco0 


Another  factor  controlling  the  fluid- structure  interaction  is  the  mass  ratio.  It  is  the  ratio 
of  the  mass  of  the  cylinder  to  the  mass  of  the  fluid  displaced  by  the  cylinder.  It  provides 
a  measure  of  the  buoyancy  effects  and  the  inertia  of  the  model  as  compared  to  the  fluid 
(Blevins  1973,  7).  It  is  given  by: 


M,» 


4  m 

npjf-s 


(39) 


Duke  University 


Page  48 


March  2006 


F49620-03- 1-0204 


Final  Report 


It  is  also  common  to  define  a  non-dimensional  vertical  displacement  coordinate. 
Therefore,  the  plunge  coordinate  can  be  divided  by  a  characteristic  length  such  as  the 
diameter  of  the  cylinder.  The  new  coordinate  is  defined  as: 

K  =  —  (40) 

D 


The  final  parameter  represents  a  ratio  of  the  fluid  frequency  to  the  structural  frequency 
and  is  given  by: 


K  = 


a>0D2 


(41) 


Using  all  these  dimensionless  parameters,  the  governing  equation  can  be  greatly 
simplified.  In  general,  the  final  equation  can  be  written  as: 


(-K2  Re2  St2  +  2  jK  Re  St  + 1  )h'x  -K  Re  C,  (Re,  St)  =  0 


(42) 


By  separating  this  equation  into  its  real  and  imaginary  parts,  a  system  of  two  equations  is 
produced.  Thomas,  et  al.  initially  demonstrated  this  method  in  2002  (Thomas,  Dowell, 
and  Hall  2002,  645).  These  equations  can  be  solved  simultaneously  for  the  two 
unknowns,  Reynolds  number  and  Strouhal  number,  for  a  specified  h\  .  The  nominal 
residual  of  the  real  and  imaginary  parts  can  then  be  written  in  the  following  form: 

'  K  2  Re2 

(-K2  Re2  St 2  +  l)/jj - —  Re^  (Re, St)] 

R(L)=  ,  2//f~  '  =0  (43) 

'  K2Re2 

(2K<^ Re St)hx - —  Im[Cz  (Re,*)] 

2/L17V  1 


where  L  is  a  vector  of  the  unknown  variables,  Re  and  St. 
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Application  of  Newton-Raphson  Technique 

The  Newton-Raphson  technique  is  then  used  to  solve  for  the  correct  combination 
of  Reynolds  number  and  Strouhal  number  for  which  the  residual  goes  to  zero.  The 
Newton-Raphson  technique  is  an  efficient  and  stable  way  to  quickly  solve  the  system  of 
equations  (Thomas,  Dowell,  and  Hall  2002,  645).  The  method  requires  the  user  to 
choose  an  initial  value  for  Reynolds  number  and  Strouhal  number.  The  method  then  uses 
the  HB  solver  to  determine  the  real  and  imaginary  parts  of  the  unsteady  lift  coefficient. 
Next,  the  Reynolds  number  is  perturbed  by  a  small  amount  and  the  change  in  the  forces 
with  respect  to  Reynolds  number  is  evaluated.  The  procedure  is  repeated  for  the  Strouhal 
number.  In  vector  notation,  this  procedure  can  be  expressed  as: 

~8R(L) 
dL 


dR  dR 


d  Re  GSt 


(44) 


where  —  and  are  given  by  the  following  relations, 
d  Re  dSt 


dR(L" )  R(L" ,  R  e"  +  s)~  R(Ln ,  Re" ) 


d  Re 


(45) 


dR(L")  '  R(Ln , St"  +  e)- R(L” , St") 
dSt  e 


(46) 


for  small  values  of  s.  The  resulting  values  can  be  put  into  matrix  fonn  and  then  the 
inverse  is  taken.  This  quantity  is  then  multiplied  by  the  original  solution  and  subtracted 
from  the  initial  guess  of  the  Strouhal  number  and  Reynolds  number.  Therefore,  this  is  an 
iterative  process  and  can  be  written  as: 

"o7?(L")n  ' 


r+1  =  Ln  - 


dL 


R(L") 


(47) 
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This  gives  a  new  update  for  Re  and  St  and  the  process  is  repeated  until  the  solution 
converges.  This  technique  is  advantageous  because  it  only  requires  a  few  iterations  to 
achieve  convergence. 

Preliminary  results  by  Jeffrey  P.  Thomas  show  that  this  is  a  viable  method  for 
predicting  the  Strouhal  number  and  Reynolds  number  of  the  self-excited  cylinder  for  a 
range  of  Reynolds  numbers.  A  plot  of  these  results  is  shown  below  in  Figure  22. 


80  90  100  110  120  130  140  150 

Reynolds  Number,  Re 


Figure  23.  Cylinder  Oscillation  Amplitude  versus  Reynolds  Number  with  Varying 
Mesh  Sizes  and  Comparison  with  the  Experimental  Data  (Thomas  2004) 

An  interesting  feature  of  this  figure  is  the  fact  that  there  are  two  different  amplitudes 

obtained  for  one  Reynolds  number,  which  is  not  possible  to  obtain  experimentally.  This 

may  be  explained  by  the  fact  that  two  solutions  are  present  -  a  stable  and  an  unstable  one. 

Therefore,  the  experiment  will  not  capture  the  unstable  solution  and  will  jump  directly  to 

the  stable  one.  In  addition,  a  plot  of  the  Strouhal  number  as  a  function  of  Reynolds 
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number  was  constructed  for  frequencies  within  the  lock-in  region.  The  results  show  that 
the  HB  method  predicts  the  relationship  between  Reynolds  number  and  Strouhal  number 
fairly  well  within  the  lock-in  region.  Therefore,  the  FIB  method  can  be  coupled  with  a 
Newton-Raphson  solver  to  accurately  model  the  aeroelastic  cylinder  case. 


100  105  110  115  120  125  130  135 

Reynolds  Number,  Re 


— HB  Method,  129x65  mesh  — ■—  Experiment 

Figure  24.  Comparison  of  Relationship  Between  Reynolds  Number  and  Strouhal 

Number  for  Aeroelastic  Case 

Experimental  Procedure  and  Results 

To  validate  these  results,  the  solution  was  compared  with  those  obtained 
experimentally  by  Anagnostopoulos  and  Bearman.  They  conducted  their  experiment  in  a 
water  channel  by  mounting  a  1 .6  mm  diameter  circular  cylinder  on  two  steel  springs  on  a 
horizontal  shaft  (Anagnostopoulos  and  Bearman  1992,  41).  Using  the  experimental 
parameter  values  given  for  the  natural  frequency,  the  damping  factor  and  the  mass  ratio 
and  calculating  K  from  given  measurements,  a  direct  comparison  can  be  made. 
Anagnostopoulos  and  Bearman  examined  Reynolds  numbers  in  and  out  of  the  lock-in 
region.  As  predicted,  they  found  that  at  a  Reynolds  number  before  lock-in,  the  vortices 
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were  shed  at  the  Strouhal  frequency  but  the  cylinder  was  oscillating  at  a  different 
frequency  so  two  frequencies  were  present  in  the  system  (Anagnostopoulos  and  Beannan 
1992,  43).  However,  when  the  Re  was  increased  to  104,  the  cylinder  oscillation 
frequency  synchronized  with  the  shedding  frequency  and  the  amplitude  of  the  cylinder 
oscillations  greatly  increased.  They  found  that  the  maximum  oscillation  amplitude 
occurred  near  the  middle  of  the  lock-in  region  (Anagnostopoulos  and  Bearman  1992,  47). 
When  the  Reynolds  number  was  increased  further  to  a  value  of  about  126,  the  shedding 
becomes  unlocked  once  again.  Therefore,  by  integrating  the  Newton-Raphson  technique 
into  the  HB  solver,  a  solution  for  the  cylinder  self-excited  aeroelastic  problem  has  been 
obtained  and  compared  with  Anagnostopoulos’  experimental  data. 

3.4  TWO-DIMENSIONAL  AIRFOIL  (Cl  CASE) 

Having  verified  that  the  technique  can  accurately  predict  the  shedding  frequency 
for  flow  about  a  cylinder,  we  next  considered  a  flow  instability  about  a  two-dimensional 
airfoil  section  of  a  modem  compressor  blade  (Cl)  operated  at  an  off-design  condition. 
First,  a  time  domain  version  of  the  flow  solver  was  generated  and  the  unsteady  normal 
force  on  the  blade  was  detennined  as  a  function  of  time  (see  Figure  25).  It  is  noted  that 
the  amplitude  has  converged  at  50  time  steps  /  iteration.  This  was  then  FFT’ed  to 
determine  the  frequency  content  with  the  results  shown  in  Figure  25  (right  side). 
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Time  Domain  Solutions 


Fourier  Transform  of  Time  Solutions 
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Figure  25.  Time  Domain  Solution  and  Fourier  Transform  of  Time  Solutions  for 
Three  Different  Timesteps/Period  (IBPA  =  0) 

To  verify  this  frequency,  the  results  were  then  compared  with  those  obtained 
using  the  HB  method.  The  blade  geometry  as  well  as  the  Mach  contours  for  the  steady 
case  is  shown  in  Figure  26.  First,  a  “steady”  flow  solution  was  obtained.  Figure  27 
shows  that  the  solution  residual  does  not  approach  zero.  Instead,  the  residual  oscillates. 
This  oscillation  is  an  indication  of  an  underlying  physical  flow  instability.  However, 
because  the  steady  flow  solver  uses  local  time  stepping  and  multi-grid  acceleration,  we 


cannot  infer  the  frequency  of  the  physical  instability  from  this  calculation. 
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Figure  26.  Mach  contours  for  the  Steady  Solver  for  2-D  Flow  About  Cl 

Compressor  Blade 


Figure  27.  Convergence  History  of  Steady  Flow  Computation  for  2-D  Flow  About 

Cl  Compressor  Blade 

Next,  the  nonlinear  harmonic  balance  solver  was  used  to  compute  the  physical 
frequency  with  interblade  phase  angle  of  zero  degrees.  The  procedure  is  to  search  for  the 
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frequency  that  allows  the  harmonic  balance  flow  solver  to  converge  to  a  zero  residual. 
Figure  28  shows  the  solution  residual  of  the  harmonic  balance  solver  as  a  function  of 
frequency  for  the  zero  interblade  phase  angle  case.  Note  that  a  frequency  of  about  1250 
Hz,  the  residual  (for  five  harmonics)  drops  dramatically  to  zero,  indicating  that  this  is  the 
frequency  of  the  flow  instability.  Note  the  very  narrow  trough  in  the  residual  curve  in 
Figure  28,  complicating  the  search  for  the  correct  frequency.  An  alternative  procedure  is 
to  search  for  the  frequency  that  produces  a  constant  phase  in  a  global  quantity  such  as  the 
first  harmonic  of  the  lift  using  the  same  method  as  described  above  for  the  cylinder. 
Figure  29  (left  side)  shows  the  phase  change  per  iteration  of  the  harmonic  balance  flow 
solver  when  the  convergence  reaches  a  steady  state.  In  this  case,  the  phase  change  is 
nearly  linear  with  frequency.  This  makes  interpolation  to  the  correct  frequency  a  simple 
matter,  and  requires  many  fewer  calculations. 
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Figure  28.  Residual  of  Harmonic  Balance  Solver  (2D  Cl  Blade,  IBPA  =  0) 

As  can  be  seen  in  Figure  29,  there  is  excellent  agreement  between  the  FFT’ed 
time  domain  results  and  the  HB  method  frequencies.  However,  it  appears  that  the  HB 
method  over  predicts  the  amplitude  of  the  first  harmonic  force  coefficient  and  under 
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predicts  the  second  harmonic.  The  higher  harmonic  amplitudes  show  good  agreement. 
Therefore,  the  HB  method  is  a  very  effective  tool  for  determining  the  NSV  frequency  but 
does  not  accurately  predict  the  amplitude  for  the  first  and  second  harmonics  for  this  2-D 
test  case. 
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Figure  29.  2D  Cl  Case:  HB  Frequency  Comparison  with  Time  Domain  Solution 

(IBP  A  =  0) 


3.5  THREE-DIMENSIONAL  CASES 


Next,  the  HB  methodology  was  applied  to  three-dimensional  real  world 
turbomachinery  blade  applications.  Based  on  input  from  industry,  NSV  was  encountered 
in  experimental  rig  testing  for  two  of  the  three  test  cases.  In  particular,  the  study 
examined  a  modern  first  stage  compressor  rotor  blade  (Cl),  a  modem  first  stage  fan  blade 
(HI),  and  a  modern  fan  vane  blade  (H2).  Although  NSV  was  not  encountered  for  the  HI 
case,  an  analysis  was  performed  to  show  that  NSV  is  not  predicted. 

Cl  CASE 


The  first  case  studied  was  a  compressor  rig  test  where  the  first  stage  modern 
compressor  blades  encountered  NSV.  For  this  test,  there  were  35  rotor  blades  and  56 
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inlet  guide  vanes  (IGV’s).  Complete  details  of  the  study  can  be  found  in  Kielb’s  paper 
entitled,  “Blade  Excitation  by  Aerodynamic  Instabilities  -  A  Compressor  Blade  Study” 
and  are  summarized  below.  The  stage  1  blades  are  observed  to  experience  a  significant 
first  torsion  mode  (IT)  response  at  lower  speed  that  shifts  to  a  second  torsion  mode  (2T) 
response  at  somewhat  higher  speed.  Figure  30  shows  a  typical  plot  of  blade  response 
frequency  and  amplitude  versus  rotational  speed.  The  vertical  lines  at  fixed  values  of 
rotor  speed  and  frequency  are  a  measure  of  blade  response  amplitude  at  that  frequency. 
The  response  at  low  speed  is  moderate  separated  flow  vibration  (SFV)  response  of  the 
first  flex  (IF)  and  IT  modes.  SFV  is  a  broadband  buffeting  response  of  the  blades.  This 
SFV  IF  response  is  followed  at  higher  speeds  by  a  significant  NSV,  which  excites  the  IT 
mode  to  a  high  level  of  response.  This  response  switches  from  a  higher  frequency  (2661 
Hz)  excitation  to  a  lower  frequency  (2600  Hz)  excitation  at  a  somewhat  higher  speed.  As 
the  speed  increases,  the  response  switches  to  a  2T  mode  excitation.  The  NSV  excitation 
of  the  IT  mode  exists  from  approximately  12700  to  12880  rpm. 
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Rotor  Speed  (rpm) 


Figure  30.  Strain  Gage  Response  of  First-Stage  Rotor  Blades  of  Compressor  Rig 

Unsteady  pressure  was  also  measured  on  the  compressor  casing  at  numerous  axial  and 
circumferential  locations.  Example  data  from  a  pressure  transducer  (located  aft  of  the 
rotor  1  blades  near  the  stage  1  vanes)  is  shown  in  Figure  31.  Significant  response  at 
frequencies  of  35 16  Hz  and  3662  Hz  was  observed  at  all  axial  and  circumferential 
measurement  locations. 
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Figure  31.  Casing  Unsteady  Pressure  Measurements 

First,  the  unsteady  CFD  code,  TURBO,  (Chen  &  Briley  (2001)),  was  used  to 
investigate  the  NSV  observed  in  the  rig  test.  A  single  row,  five  passage,  mesh  modeled 
one-seventh  of  the  rotor  circumference.  The  mesh,  consisting  of  188  axial,  56  radial,  and 
280  circumferential  grid  points,  contained  approximately  three  million  grid  points.  Using 
TURBO,  the  analysis  took  about  six  months  to  run  on  a  single  processor.  The  local 
unsteady  static  pressures  from  the  periodic  CFD  solution  were  investigated  at  many 
locations  on  the  blade  surface  near  the  blade  tip.  This  is  shown  in  Figure  32,  where  the 
unsteady  static  pressures  in  passage  three,  near  the  leading  edge  on  the  pressure  side,  are 
presented.  As  can  be  seen,  the  unsteady  pressure  content  at  2365  Hz  and  4370  Hz  is 
much  greater  than  that  at  the  vane  passing  frequency.  The  predicted  NSV  frequency 
(2365  Hz)  is  approximately  9%  lower  than  that  measured  in  the  rig  test.  The  TURBO 
results  show  that  NSV  are  primarily  a  coupled,  suction  side  vortex  unsteadiness  (near 
75%  span)  and  tip  flow  instability. 
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Figure  32.  Unsteady  Static  Pressure  Near  Tip  Pressure  Side  Near  Leading  Edge 

In  an  attempt  to  replicate  these  results,  similar  conditions  were  used  and  the 
analysis  was  performed  using  Duke  University’s  time  domain  and  HB  computational 
methods.  The  blade  was  excited  at  approximately  2.6  Khz.  A  time  domain  CFD 
analyses  was  conducted  using  the  running  clearance  of  0.020  in.,  zero  interblade  phase 
angle,  and  pexit/pexitCAFD  =  1.000  where  CAFD  is  the  exit  pressure  for  which  NSV  was 
encountered  in  the  experiment.  An  NSV  frequency  of  2.9  Khz  was  predicted  (see  Figure 
33).  We  next  used  our  harmonic  balance  method  on  this  case.  HB  solution  stability 
problems  were  encountered  when  running  at  the  NSV  speed.  The  phase  change  per 
iteration  of  the  unsteady  torque  is  shown  in  Figure  34  for  five  assumed  frequencies  with 
the  same  conditions  as  the  time  domain  case.  As  can  be  seen  the  phase  change  reaches  a 
constant  value  for  2.5  and  2.6  KHz.  For  2.6  and  2.75  KHz  the  phase  change  does  not 
converge  to  a  constant  value.  This  is  likely  because  the  fluid  dynamic  instability  contains 
multiple  (irrational)  frequencies. 
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Figure  33.  3D  Cl  Case:  FFT  of  Time  Domain  Solution  (IBPA  =  0,  running 

clearance,  pexit/pexitCAFD  =  1.000) 


Figure  34.  3D  Cl  Case:  Phase  Change  Per  Iteration  (IBPA  =  0,  0.020  in.  Running 

Tip  Clearance,  pexit/pexit  ,exp  1.000) 

For  this  study,  a  0.020  inch  tip  clearance  was  used  in  accordance  with  the 
experimental  data  and  computational  data  obtained  from  TURBO.  Due  to  the  limitations 
of  the  HB  method  as  discussed  above,  only  the  time  domain  results  are  presented  for  the 
running  clearance  case.  The  two  main  design  parameters  that  were  varied  were  the  back 
pressure  and  the  tip  clearance.  For  the  0.020  inch  case,  two  frequencies  are  present. 
Flowever,  one  of  the  frequencies  happens  to  be  very  low.  In  Figures  35  and  36,  an  FFT 
of  the  time  domain  results  for  the  blade  torque  is  presented  for  back  pressures  of  1.015 
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and  1.020.  Figure  36  demonstrates  the  results  in  the  0-100  Hz  range,  where  a  very  low 
frequency  component  is  present  in  the  solution.  Furthermore,  the  time  domain  results  for 
the  mass  flow  rates  for  back  pressures  of  1.010,  1.015,  1.020  are  presented  in  Figure  37 
and  compared  with  the  results  obtained  experimentally  (CAFD)  and  using  TURBO.  It  is 
noted  that  the  flow  is  steady  for  pexit/pexitCAFD  =  1.010.  Finally,  Figure  38  shows  the  time 
history  results  for  various  back  pressures  for  the  running  tip  clearance  case. 
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Figure  35.  FFT  of  Blade  Torque  for  Two  Different  Back  Pressures  (0.020  in. 
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Figure  36.  FFT  of  Blade  Torque  for  Two  Different  Back  Pressures  (0.020  in. 
Running  Tip  Clearance)  in  Low  Frequency  Region 
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Figure  37.  Time  Domain  Solution  Mass  Flow  Rates  for  Different  Back  Pressures 

(0.020  in.  Running  Tip  Clearance) 
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Figure  38.  Time  History  of  Blade  Torque  for  Various  Back  Pressures  (0.020  in. 

Running  Tip  Clearance) 

Due  to  the  limitations  of  the  HB  method  to  only  handle  one  irrational  frequency,  a 
tight  tip  clearance  case  was  also  examined.  Figure  39  shows  time  histories  of  blade 
torque  for  the  tight  tip  clearance  for  two  different  back  pressures  and  Figure  40  shows  the 
FFT's  of  the  time-domain  solutions  along  with  a  HB  result.  As  can  be  seen,  HB  shows 
good  agreement  for  the  lower  back  pressure  case  since  there  appears  to  be  only  a  single 
frequency,  and  higher  harmonics  thereof,  present  in  the  time  domain  solution.  However, 
the  higher  back  pressure  case  has  two  distinct,  irrational  frequencies  (beats)  present,  and 
HB  can  only  handle  one  fundamental  frequency  and  its  harmonics.  A  summary  of  the 
NSV  frequency  results  obtained  can  be  found  in  Table  2.  Some  possible  solutions  to 
overcome  the  HB  method’s  inability  to  handle  more  than  one  fundamental  frequency  are 
to  extend  the  HB  method  to  handle  multiple  irrational  frequencies,  use  enforced  motion 
to  cause  “lock-in”  to  the  blade  frequency,  or  use  the  aeroelastic  solution  to  cause  “lock- 
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in”  to  the  non-linear  flutter  frequency  (very  near  the  blade  frequency).  All  of  these 
options  are  the  subject  of  current  research  study. 
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Figure  39.  Time  History  of  the  Blade  Torque  for  Two  Different  Back  Pressures 

(0.002  in.  Tight  Tip  Clearance) 
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HI  CASE 

Next,  a  study  was  conducted  for  a  first  stage  modern  fan  blade  (known  as  HI)  that 
encountered  flutter  and  not  NSV.  The  GUIde  Consortium  members  asked  that  an 
analysis  be  conducted  to  show  that  NSV  is  not  predicted.  As  a  result,  both  Euler  and  NS 
steady  and  unsteady  analyses  were  performed.  The  steady  analysis  showed  no  signs  of 
NSV.  A  linear  flutter  analysis  was  then  conducted  with  the  results  compared  to  time 
domain  results  (TURBO)  and  is  shown  in  Figure  41.  There  is  reasonable  agreement 
between  the  time  and  frequency  domain  results.  The  latest  efforts  involved  reducing  the 
back  pressure  to  produce  negative  aerodynamic  damping. 


Duke  University 


Page  67 


March  2006 


F49620-03- 1-0204 


Final  Report 


HI  Damping  VS.  ND 


Figure  41.  Time  Domain  and  HB  Results  for  Damping  as  a  Function  of  Nodal 

Diameter  for  the  HI  Case 

H2 CASE 

The  final  case  that  was  studied  was  a  modern  fan  vane  blade  known  as  the  H2 
case.  This  blade  encountered  NSV  in  experimental  rig  testing.  An  analysis  was 
performed  with  TURBO  and  it  showed  good  agreement  with  the  experimental  data. 

Since  there  is  no  clearance  on  a  vane  blade,  this  was  thought  to  be  a  good  test  case  for  the 
HB  method.  However,  the  same  problems  were  encountered  as  in  the  C 1  case.  Multiple 
irrational  frequencies  were  present  and  the  current  HB  method  can  only  handle  one 
fundamental  frequency  and  its  harmonics.  Therefore,  this  case  is  still  under 
investigation. 

4.  SUMMARY  AND  CONCLUSIONS 

In  recent  years,  new  aeromechanic al  problems  have  been  encountered  in 
turbomachinery.  In  particular,  non-synchronous  vibrations  in  turbine  blades  have  been 
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observed  by  many  engine  companies  and  occur  as  a  result  of  a  fluid  dynamic  instability 
such  as  shedding.  If  the  oscillations  become  large  enough,  blade  fatigue  and  failure  can 
occur.  It  is  difficult  to  predict  because  there  is  little  knowledge  of  the  phenomenon.  As  a 
result,  there  is  a  strong  motivation  to  develop  a  fast  computational  method  to  be  able  to 
understand  this  behavior  in  the  design  stage. 

For  this  study,  many  test  cases  were  examined.  As  an  initial  demonstration  of  the 
method,  the  flow  over  a  cylinder  was  investigated.  This  serves  as  a  useful  test  case  for 
modeling  NSV  because  it  is  a  well-studied  phenomenon  and  a  significant  amount  of 
experimental  data  is  available.  Next,  this  technique  was  extended  to  study  the  flow  over 
a  two-dimensional  airfoil.  In  particular,  we  studied  the  2D  Cl  case.  Finally,  the 
hannonic  balance  method  was  extended  to  three-dimensional  cases  from  industry  and 
compared  with  existing  time  domain  solutions  and  experimental  data.  In  particular,  it 
was  applied  to  the  Cl  case  (a  front  compressor  blade),  the  HI  case  (a  fan  blade),  and  the 
H2  case  (a  fan  vane).  Therefore,  there  are  numerous  test  cases  available  to  test  the 
hannonic  balance  method  and  establish  it  as  a  valuable  design  tool  for  future  use  by 
engine  manufacturers. 

The  numerical  approach  involved  the  use  of  a  two-dimensional  nonlinear  Navier- 
Stokes  unsteady  harmonic  balance  method  developed  by  Hall,  et  al.  One  of  the  major 
advantages  of  this  method  is  that  it  requires  one  to  two  orders  of  magnitude  less 
computational  time  than  conventional  time  marching  CFD  techniques.  In  addition,  a 
unique  phase  error  method  was  developed  to  detennine  the  precise  NSV  frequency  in 
only  a  few  hannonic  solution  computations. 
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Next,  the  HB  method  was  applied  to  the  2D  Cl  case  from  industry  as  well  as  the 
3D  Cl,  HI,  and  H2  cases.  Preliminary  results  showed  good  agreement  for  the  2D  case 
and  promising  results  from  the  3D  cases.  For  the  3D  Cl  case,  two  main  design 
parameters  were  varied  -  the  tip  clearance  and  the  back  pressure.  Using  the  running  tip 
clearance  of  0.020  in.,  two  irrational  frequencies  are  present  in  the  solution  so  the  current 
HB  method  could  not  be  used.  However,  time  domain  simulations  were  performed  to 
determine  the  NSV  frequency  for  three  different  back  pressures.  In  an  attempt  to  justify 
the  merits  of  the  HB  method,  a  tighter  tip  clearance  of  0.002  in.  was  studied.  The  HB 
method  showed  good  agreement  with  the  time  domain  results  for  the  case  of  PR  =  1.000 
because  only  a  single  frequency  and  its  hannonics  were  present  for  this  case.  Therefore, 
current  research  is  focused  on  the  following  three  methods  to  solve  for  the  NSV 
frequency.  One  solution  is  to  extend  the  HB  method  to  handle  multiple  irrational 
frequencies.  Other  solutions  are  to  use  enforced  motion  to  cause  “lock-in”  to  the  blade 
frequency  or  the  aeroelastic  solution  to  cause  “lock-in”  to  the  non-linear  flutter  frequency 
(very  near  the  blade  frequency)  as  was  done  in  the  study  of  the  cylinder. 

In  addition,  the  HI  and  H2  cases  were  studied.  For  the  HI  case,  the  steady 
analysis  showed  no  signs  of  NSV.  A  linear  flutter  analysis  was  then  conducted  and  the 
results  were  compared  to  time  domain  results  (TURBO).  Reasonable  agreement  was 
achieved  between  the  time  and  frequency  domain  results.  The  latest  research  efforts  have 
involved  reducing  the  back  pressure  to  produce  negative  aerodynamic  damping.  Similar 
to  the  Cl  case,  the  H2  case  encountered  multiple  irrational  frequencies  present  in  the 
solution.  Since  the  current  HB  method  can  only  handle  one  fundamental  frequency  and 
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its  hannonics,  the  methods  described  above  are  being  implemented.  This  case  is  still 
under  investigation. 

Therefore,  the  harmonic  balance  method  has  been  established  as  a  viable  tool  to 
solve  highly  nonlinear  unsteady  shedding  flow  over  turbine  blades,  the  ultimate  goal  of 
the  research  effort.  The  results  are  consistent  with  those  obtained  by  TURBO,  a  high 
fidelity  computational  tool.  In  addition,  the  HB  method  is  able  to  predict  both  the  NSV 
frequency  and  amplitude  for  the  cylinder  case  and  the  2-D  Cl  airfoil  case  but  only  for 
some  3-D  flows.  This  is  due  to  the  presence  of  multiple  NSV  frequencies,  which  the 
current  HB  method  cannot  handle.  Research  is  underway  to  address  this  problem  and 
possible  solution  techniques  are  presented  above.  This  method  is  a  valuable  design  tool 
because  the  frequency  can  be  added  to  the  Campbell  diagram  to  identify  regions  where 
the  vortex  shedding  frequency  may  interact  with  the  blade  natural  frequencies. 
Furthermore,  it  will  allow  engineers  to  better  understand  NSV  behavior  and  to  predict  its 
occurrence  in  the  design  stage  for  flow  over  turbine  engine  blades. 
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